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EXECUTIVE  SUMMARY 


This  report  describes  the  theory  support  of  DTRA’s  Plasma  Radiation  Source  (PRS)  program 
carried  out  by  NRL’s  Radiation  Hydrodynamics  Branch  (Code  6720)  in  FY  2000.  Included  is  work 
called  for  in  DTRA  MIPR  00-2022  -  “Plasma  Radiation  Theory  Support”  and  in  Sandia  National 
Laboratories  document  number  8591  -  “Diagnostic  and  Spectroscopic  Calculations  and  Analysis  . 

This  was  a  year  for  testing  K-shell  yield  scaling  with  atomic  number  on  the  Z  machine  at 
Sandia  National  Laboratories.  Z’s  capability  of  driving  a  PRS  load  with  15-20  mega-amperes  of 
current  allowed  us  to  explore  higher  Z  plasma  radiation  source  (PRS)  loads  like  nickel  (z=28)  and 
copper  (z=29).  The  experimental  nickel,  and  copper  K-shell  yields  of  16  kJ  and  12  kJ,  respectively, 
are  world  records.  Another  world  record  yield  was  also  set  on  Z  for  argon,  275  kJ.  This  was  a  major 
accomplishment  considering  that  it  was  the  first  time  a  gas  puff  experiment  had  ever  been  fielded 
on  the  Z  machine.  This  record  argon  yield  was  large  enough  to  be  in  the  so  called  “efficient”  scaling 
regime  and  it  represents  the  first  time  an  efficient  argon  yield  was  ever  achieved  using  a  large  radius 
nozzle.  Given  this  success,  it  would  be  interesting  to  further  explore  the  atomic  number  scaling 
of  K-shell  emission  by  imploding  a  krypton  gas  puff  load  on  Z.  These  recent  record  yields  along 
with  previous  titanium  and  stainless  steel  yields  give  us  confidence  that  the  basic  ideas  about  PRS 
behavior  that  have  been  developed  over  the  years  by  Code  6720  are  valid  and  provide  DTRA  with 
a  reasonable  foundation  upon  which  future  PRS  experiments  can  be  designed  and  built. 

Given  that  it  may  be  some  time  before  the  Z  machine  is  upgraded,  and  given  that  the 
Decade  Quad  will  be  only  a  5  MA  to  8  MA  machine,  it  will  be  awhile  before  K-shell  yields 
are  substantially  increased  or  substantially  extended  in  energy  and  efficiency  due  to  higher  driving 
currents.  Therefore,  in  order  to  increase  yields  with  our  present  inventory  of  machines  supplying 
the  drive  current,  it  becomes  increasingly  important  to  fine  tune  our  PRS  knowledge  so  that  we 
can  design  higher  performance  PRS  loads.  This  fine  tuning  has  been  initiated  by  substantial 
progress  in  many  aspects  of  PRS  technology.  In  this  report  the  progress  made  in  load  design, 
model  development  and  assessment,  and  in  diagnostic  development  for  the  purpose  of  fine  tuning 
our  understanding  of  Z-pinch  behavior  is  discussed  section  by  section. 

(1)  Experience  has  taught  us  that  improved  symmetry  and  stability  are  the  keys  to  better 
load  design.  Indeed,  it  was  close  consideration  of  these  key  issues  that  allowed  Sandia  National 
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Laboratories  to  achieve  their  record  radiative  powers  with  wire  array  loads  on  the  Saturn  machine, 
which  ultimately  led  to  dramatic  changes  in  their  entire  fusion  program.  A  proven  method 
for  achieving  better  initial  load  symmetry  is  to  use  more  wires.  However,  the  drive  power  of 
the  machine  sets  limits  on  the  amount  of  mass  that  can  be  imploded  and  ionized  to  the  K- 
shell.  This  in  turn  limits  the  number  of  wires  that  can  be  used  in  an  experiment  because  of  the 
practical  limitations  on  the  minimum  mass-per-unit-length  of  wire  that  can  be  manufactured.  In 
order  to  (1)  improve  comparisons  between  theory  and  experiment,  (2)  avoid  the  problems  and 
uncertainties  associated  with  wire  gap  spacing,  or  for  gas  puffs,  uncertainties  in  the  initial  gas 
distribution,  (3)  test  our  understanding  of  radiation-hydrodynamics,  (4)  explore  alternate  avenues 
for  initiating  plasmas  and  (5)  determine  prospects  for  producing  high  density  compressions  and  K- 
shell  blackbody  radiators,  we  are  proposing  a  new  family  of  loads  that  are  thin  cylinders  of  plastic 
foil.  These  cylindrical  foils  allow  for  the  possibility  of  depositing  thin  layers  of  higher  atomic 
number  materials  onto  the  foil  surface  that  become  the  K-shell  emitters  following  implosion.  In 
a  sense  this  proposed  load  represents  an  evolutionary  step  along  the  way  toward  x-ray  simulators 
with  enough  current  to  implode  solid  cylinders  with  sufficient  amounts  of  high  Z  materials  to 
achieve  efficient  K-shell  emission.  This  load  design  is  discussed  in  detail  in  Section  I. 

(2)  There  has  to  be  reliable  feedback  between  experiment,  model  development  and  load  design 
if  there  is  to  be  sustained  progress  in  understanding  Z-pinch  behavior.  Plasma  diagnostics  and 
analysis  are  responsible  for  a  large  part  of  this  feedback.  By  directly  comparing  as  much  spectral 
information  as  is  available  (e.g,  measured  total  powers,  K-shell  powers,  the  size  of  the  K-shell 
emission  region,  K-shell  line  ratios,  spatially  resolved  spectra,  and  time  integrated  and  resolved 
spectra)  with  corresponding  calculated  information  using  our  collisional  radiative  equilibrium 
ionization  dynamics  models  and  our  one-dimensional  (ID)  radiation  transport  plasma  models,  we 
are  making  steady  progress  in  improving  our  inferences  about  the  state  of  the  plasma.  Our  models 
have  been  improved  to  the  point  where  time  integrated  and  time  resolved  radial  gradients  in  plasma 
temperature  and  density  can  be  inferred.  The  work  that  is  presented  in  Section  II  of  this  report 
demonstrates  how  these  models  can  be  used  to  infer  the  radial  profiles  of  electron  temperature  and 
density  by  matching  model  calculations  with  radially  resolved  spectroscopic  line  measurements 
of  the  Z  pinches  created  in  Double  Eagle  aluminum-silicon  wire  array  experiments.  Diagnostic 
techniques  and  analysis  such  as  this  provide  us  with  knowledge  of  the  state  of  the  experimental 
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plasmas.  Matching  the  state  of  the  plasma  as  calculated  using  our  ID  and  2D  radiation  MHD 
models  to  this  inferred  knowledge  of  the  experimental  state  of  the  plasma  is  the  main  method  we 
employ  to  improve  our  understanding  of  the  dynamics  of  PRS  loads. 

(3)  A  straightforward  approach  to  establishing  a  theoretical  relationship  between  the  initial 
asymmetries  or  the  Rayleigh  Taylor  (RT)  instabilities  of  a  z-pinch  and  its  ability  to  efficiently 
radiate  K-shell  photons  would  be  to  perform  two-dimensional  (2D)  or  three-dimensional  (3D),  time 
dependent,  non-local  thermodynamic  equilibrium  (non-LTE),  full  radiation  transport,  magneto¬ 
hydrodynamic  (MHD)  calculations.  However,  this  capability  has  never  been  demonstrated  by  the 
z-pinch  community  because  the  computer  requirements  for  performing  the  zone-to-zone  coupling 
needed  to  solve  the  equation  of  radiative  transfer  in  a  2D  or  a  3D  plasma  are  severe.  In  the 
past,  our  2D  MHD  model,  PRISM,  which  stands  for  Plasma  Radiating  Imploding  Source  Model, 
inadequately  modeled  the  radiation  field  by  assuming  either  that  the  plasma  was  optically  thm 
or  that  a  local  escape  model  could  be  used  to  model  some  of  the  photon  self-absorption  physics. 
Another  model  widely  used  by  the  z-pinch  community,  the  LTE  radiative  diffusion  model,  is  also 
inadequate  for  modeling  the  radiation  physics  of  PRS  loads  designed  for  K-shell  emission  because 
neither  the  diffusion  approximation  nor  the  LTE  assumption  are  valid  for  K-shell  plasma  conditions 
produced  by  present  day  Z-pinches.  In  Section  HI  we  describe  a  newly  developed,  computationally 
efficient,  multi-dimensional,  radiation  transport  method  for  modeling  Z-pinch  plasmas  that  are  in 
collisional  radiative  equilibrium.  This  method  is  superior  to  diffusion  methods  because.  (1)  the 
source  function  is  based  on  a  self-consistent  non-LTE  calculation,  not  an  ad  hoc  assumption  that  it 
is  Planckian,  (2)  the  optically  thin  radiation  is  allowed  to  escape,  not  subjected  to  an  inappropriate 
diffusion  equation,  and  (3)  accurate  line-by-line  escape  probabilities  are  used,  not  diluted  as  in 
diffusion,  by  frequency  groups  that  mix  line  and  continuum  opacities.  This  multi-dimensional 
radiation  transport  model  is  employed  by  our  2D  PRISM  code  to  analyze  the  effect  of  Rayleigh- 
Taylor  instability  growth  on  argon  K-shell  yield  production.  In  particular  gas  puff  loads  driven  by 
a  Z  machine  circuit  are  modeled. 

(4)  Wire  arrays  have  proven  to  be  reliable  PRS  loads  for  achieving  record  radiative  powers 
and  yields.  The  progress  made  in  the  research  and  development  of  these  loads,  which  represent 
such  a  prominent  element  of  DTRA’s  PRS  program,  is  described  in  Section  IV.  Early  insight  into 
the  behavior  of  a  wire  array  load  was  provided  by  the  Wire  Dynamic  Model  (WDM)  wherein 
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each  wire  is  treated  as  a  distinct  circuit  element.  This  model  has  successfully  described  mass  and 
current  interchange  in  nested  arrays  over  broad  current  and  mass  ranges.  Improvements  to  the 
WDM  model  that  allow  the  treatment  of  bias  inductance  and  wire  resistance  are  discussed  as  well 
as  treatment  of  the  precursor  plasma  that  forms  as  the  wires  heat.  Treating  this  tenuous  precursor 
plasma  stresses  conventional  MHD  codes  to  the  point  of  validity.  In  order  to  more  accurately 
model  the  precursor  plasma  we  are  evolving  the  WDM  model  into  a  true  fluid  particle  model  and 
progress  made  toward  this  goal  is  described. 

(5-6)  Extensive  experience  in  modeling  plasma  radiation  sources  over  many  years  has  raised 
our  awareness  of  the  large  sensitivity  of  model  results  to  initial  conditions,  zoning  assumptions, 
model  platforms  -  i.e.,  Eulerian  versus  Lagrangian  reference  frames,  and  transport  coefficients.  In 
the  work  presented  in  Section  V,  sensitivity  issues  related  to  the  design  and  modeling  of  titanium, 
nested-plasma,  Z-pinch  loads  are  investigated.  Specifically,  the  K-shell  emission  sensitivity  to  the 
initial  plasma  mass  distribution,  and  the  role  of  electron  viscosity  in  directly  converting  viscous 
energy  into  electron  thermal  energy  are  examined.  Typically,  in  past  calculations,  the  electron 
viscosity  was  ignored  and  all  the  the  viscous  energy  went  into  heating  the  ions.  Consequently, 
the  calculated  ion  temperatures  are  very  high  -  higher  than  line  broadening  measurements  have 
indicated.  In  order  to  examine  the  sensitivity  of  K-shell  emission  to  this  viscosity  issue,  the 
standard  form  of  artificial  viscosity  was  used,  but  the  heating  it  produced  was  divided  between 
electrons  and  ions'according  to  the  ratio  of  their  classical  viscosities.  The  role  of  electron  viscosity 
in  plasma  shock  heating  is  discussed  in  detail  in  Section  VI. 

(7)  Decontamination  of  the  anode  in  the  bremsstrahlung  diode  has  substantially  improved 
bremsstrahlung  radiation  production.  This  improvement  was  demonstrated  experimentally  on 
Gamble  II  last  year  and  a  similar  set  of  experiments  is  to  take  place  on  DECADE  MODULE  2 
(DM2)  in  March  of  2001.  The  success  of  the  DM2  decontamination  experiments  is  important  for 
achieving  similar  improvements  in  bremsstrahlung  production  on  the  DECADE  machine.  This 
section  describes  experimental  and  theoretical  work  that  has  been  and  will  be  done  to  understand 
the  role  that  heating  the  tantalum  anode  in  a  bremsstrahlung  diode  to  2000  -  3000  K  plays  in 
improving  bremsstrahlung  production. 


I.  A  LOAD  DESIGN  PROPOSAL 


Three  important  features  of  Z-pinch  dynamics  that  influence  its  PRS  performance  are  evident 
from  past  comparisons  between  Z-pinch  theory  and  Z-pinch  experiments.  One,  experiments  using 
tracer  elements  as  well  as  theory  show  that  K-shell  emissions  originate  from  the  core  of  the  plasma 
that  assembles  on  axis,  and  this  radiation  generally  traverses  a  colder  (L-shell  emitting)  outer 
plasma.  Two,  the  outer  portion  of  the  assembled  plasma  is  generally  much  larger  than  the  core 
region,  and  it  extends  well  beyond  one-dimensional  (1-D)  computer  calculation  predictions.  Part 
or  possibly  all  of  this  added  extent  has'  sometimes  been  explained  by  the  growth  of  Rayleigh- 
Taylor  instabilities  at  the  outer  surface  of  the  plasma.  However,  this  is  not  a  firm  comparison 
between  theory  and  experiment  since  it  depends  on  how  much  initial  perturbation  one  applies  at  the 
beginning  of  a  calculation  and  on  how  much  one  imagines  that  the  lengthy  spikes  being  calculated 
correspond  to  the  qualitative  features  that  are  seen  in  x-ray  pinhole  pictures.  In  any  case,  both  1-D 
and  2-D  (r-z)  calculations  assume  perfectly  cylindrically  shaped  plasmas  initially  that  do  not  exist 
in  experiments.  Consequently,  three,  experiments  always  show  a  dependence  of  K-shell  yields  and 
powers  on  initial  wire  gap  spacing  or  initial  gas  distribution.  The  difficulty  of  controlling  the  initial 
mass  distribution  of  the  plasma  load  generally  degrades  the  ability,  in  particular,  of  moderate  atomic 
number  plasmas  to  compress  into  tight  uniform  pinches  on  axis  that  can  radiate  with  reasonable 
efficiency  from  the  K-shell. 

The  constraint  that  has  prevented  elimination  of  this  unwanted  dependence  on  wire  number 
and  wire  gap  spacing  in  the  past  has  always  been  the  limited  amount  of  current  that  can  be  delivered 
to  a  Z-pinch  load.  With  too  low  a  current  drive,  the  masses  of  cylindrical,  foil-type  loads  made 
from  solid  materials  are  generally  too  large  to  be  imploded  to  high  enough  velocities  to  radiate  well, 
if  at  all,  from  the  K-shell.  However,  the  Z  generator  at  Sandia  National  Laboratories  has  reduced 
this  limitation  considerably  and  may  allow  for  the  fabrication  of  a  new  family  of  loads  of  the  kind 
proposed  here.  Figure  1  illustrates  a  possible  design.  A  set  of  support  pins  are  mounted  at  the  initial 
radius  of  the  load  on  both  the  nominal  cathode  and  anode  surfaces.  For  simplicity,  the  figure  shows 
only  six  pins  and  the  same  number  of  support  pins  as  return  current  posts.  They  are  situated  at  the 
same  angular  positions  as  the  posts  to  compensate  somewhat  for  the  only  approximately  cylindrical 
shape  of  the  surface.  A  thin  foil  of  plastic  is  stretched  and  attached  to  the  support  pins  and  becomes 
the  Z-pinch  load  that  undergoes  implosion.  As  shown  in  the  figure,  this  load  design  allows  for  the 
possibility  of  depositing  thin  layers  of  higher  atomic  number  materials,  such  as  titanium,  iron,  or 
copper,  onto  the  inner  surface  of  the  foil  that  become  the  K-shell  emitters  following  implosion. 

A  Z-pinch  load  of  the  kind  depicted  in  Figure  1  solves  a  number  of  problems  and  creates 
others.  The  use  of  support  pins  gives  the  load  structural  stability  unavailable  in  an  unsupportable, 
thin  film,  design  (although  foam  cylinders  may  have  such  stability  and  sufficiently  light  weight 
when  they  become  available).  When  Z-pinch  experiments  were  done  in  the  past  with  thick  wires,  it 
was  observed  that  some  fraction  of  the  wires  had  been  left  behind  unimploded  and  in  their  original 
positions  at  the  time  when  the  pinch  was  assembling  on  axis.  Thus,  the  support  pins  should  not 
interfere  with  the  foil  implosion  if  they  are  not  allowed  to  form  a  closed  electrical  circuit.  The  use 
of  a  foil  eliminates  gap  spacings;  however,  there  are  problems  (as  with  wires)  as  to  how  one  should 
initiate  a  discharge  in  such  a  way  as  to  convert  solid  material  into  a  cylindrically  symmetric  plasma. 
However,  if  such  initial  conditions  can  be  achieved,  the  first  truly  unambiguous  comparisons 
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between  1-  and  2-D  calculations  and  experiment  should  be  possible  without  the  need  to  guess  at  the 
influence  of  hidden  experimental  variables.  By  eliminating  gap  spacings,  one  could  also  eliminate 
the  premature  arrival  of  plasma  on  axis  ahead  of  the  main  body  of  plasma.  Calculations  show 
that  such  precursor  plasma  is  preheated,  and  that  it  impedes  the  compression  of  the  plasma  core, 
which  is  beneficial  to  the  production  of  K-shell  radiation.  The  implosion  of  a  plastic  foil  could 
become  the  first  implosion  of  a  true  plasma  shell.  It  would  allow  a  far  more  accurate  evaluation 
of  the  influence  of  Rayleigh-Taylor  instabilities  on  the  K-shell  performance  of  Z-pinch  loads  than 
has  been  possible  to  date.  Moreover,  the  fact  that  the  load  is  not  initially  in  the  shape  of  a  perfect 
cylinder  could  actually  be  a  benefit  since  it  compensates  somewhat  for  the  slight  reduction  in  the 
size  of  the  j  x  B  forces  away  from  the  support  rods  because  of  the  finite  number  of  return  current 
posts  generally  used  in  experiments. 

The  ability  to  undercoat  the  plastic  with  angstrom-sized  thicknesses  of  moderate  atomic 
number  ( Z )  elements  could  also  be  a  major  benefit.  Calculations  have  shown  (see  Section  xx  of 
this  final  report),  that  in  a  true  1-D-like  implosion,  the  plasma  core  is  compressed  and  heated,  while 
the  compressing  plasma  remains  at  higher  densities  and  lower  temperatures.  The  K-shell  emission 
in  moderate  Z  element  implosions  comes  from  this  core  plasma.  However,  it  is  not  necessary 
to  have  a  moderate  Z  element  do  the  compressing  and  heating  of  the  core.  The  replacement  of 
moderate- Z  with  a  low-Z  element  in  the  outer  plasma  has  several  advantages.  Plastic  is  less  dense 
than  titanium  or  copper,  and  therefore,  a  thicker  shell  of  material  can  be  used  in  the  load  design. 
The  density  of  plastic  is  roughly  1  gm/cm3,  4.5  times  less  dense  than  titanium  and  9  times  less 
dense  than  copper.  Thus,  a  1  cm  radius  cylinder  of  10  /xm  thick  plastic  would  only  have  a  mass 
of  6.3  mg/cm,  while  an  identical  titanium  cylinder  would  be  4.5  times  heavier  and  an  identical 
copper  cylinder  9  times  heavier.  Furthermore,  by  having  a  fully  ionized  carbon-hydrogen  plasma 
in  front  of  the  K-shell  emitting  core  plasma  instead  of  a  partially  ionized,  moderate- Z,  plasma,  the 
absorption  losses  of  the  K-shell  emissions  are  reduced. 

Two  sets  of  0-D  calculations  were  carried  to  illustrate  the  potential  interest  of  experiments  on 
the  Z  generator  that  use  pure  plastic  loads  like  the  one  in  Figure  1.  Figures  (2a)  and  (2b)  contain 
the  results  of  these  calculations.  For  reference  purposes,  Figure  (2a)  shows  the  locations  of  the 
mass  breakpoint  curves  (both  soft  and  hard)  for  a  carbon  plasma  as  well  as  the  implosion  velocities 
corresponding  to  r\  values  of  7/  =  4  and  77  =  20  for  carbon.  For  implosion  velocities  between  or 
greater  than  these  values,  carbon  is  a  highly  optically-thick  and  efficient  emitter  of  K-shell  radiation. 
Two  series  of  0-D  implosion  calculations  of  different  massed  carbon  loads  were  carried  out  using 
an  early  (somewhat  approximate  but  close)  version  of  the  Z  generator  circuit.  In  one  series  of 
calculations,  the  initial  load  radius  was  1  cm,  and  in  the  other  series,  it  was  1.5  cm.  Figure  (2b) 
shows  the  implosion  velocities  these  loads  achieved  in  calculations  that  were  stopped  after  the  load 
had  imploded  to  a  radius  of  2  mm.  They  demonstrate  that  the  Z  generator  has  the  capability  to 
drive  carbon  loads  of  masses  between  10  and  20  mg/cm  to  velocities  well  beyond  those  needed  for 
the  efficient  production  of  (highly  opaque)  K-shell  radiation.  Important  questions  thus  arise  that 
are  in  need  of  experimental  investigation.  How  close  to  a  blackbody  spectrum  can  such  emissions 
get?  How  close  to  the  tungsten  emission  spectra  could  the  carbon  K-shell  spectra  get  as  a  result  of 
the  greater  mass  of  the  carbon  loads  and  their  potentially  greater  implosion  symmetry? 

A  previous  theoretical  study1  of  CH  foams  heated  by  ion  beams  from  the  Particle  Beam 
1  J.  P.  Apruzese,  J.  E.  Rogerson,  R.  W.  Clark,  and  M.  C.  Coulter,  Naval  Research  Laboratory 
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Fusion  Accelerator  II  at  Sandia  National  Laboratories  (Z’s  predecessor)  addressed  some  of  these 
questions.  In  this  study,  it  was  demonstrated  that  low  atomic  number  plasmas  of  the  mass  range 
and  material  composition  discussed  above  can  radiate  from  the  K-shell  as  near-blackbody  sources 
(see  Figures  3  through  7).  Two  factors  contribute  to  this  behavior.  First,  the  extreme  optical 
depths  produced  by  such  large  mass  plasmas  saturates  the  line  emission  at  the  blackbody  curve 
and  increases  continuum  emissions.  Second,  low  atomic  number  elements  such  as  carbon  have 
high  collisional  excitation  and  deexcitation  rates  coupling  different  ionization  stages  and  bound 
states  within  a  given  stage.  This  coupling  results  in  local  thermodynamic  equilibrium  (LTE)  being 
achieved  at  lower  ion  densities  than  is  the  case  for  medium  or  high  atomic  number  elements. 

Figures  3-7  show  the  progression  towards  blackbody  emission  of  different  amounts  of  carbon 
plasma  at  the  uniform  plasma  conditions  indicated  in  the  figures.  Two  sets  of  calculations  were 
done.  In  one  set,  the  plasma  had  a  uniform  temperature  of  100  eV;  in  the  other,  its  temperature  was 
raised  to  200  eV.  In  all  cases,  the  plasma  was  a  cylinder  of  4  mm  diameter.  At  100  eV,  5  mg/cm  of 
carbon  plasma  was  already  sufficient  to  closely  approximate  blackbody  emission  [Fig.  (3)],  while 
10  mg/cm  fully  achieved  blackbody  emission  [Fig.  (4)].  For  a  cylinder  2  cm  in  length  with  a  4 
mm  diameter  radiating  as  a  100  eV  blackbody,  the  power  output  would  be  26  TW,  a  value  easily 
achievable  on  the  Z  generator.  By  upping  the  blackbody  temperature  to  200  eV,  the  power  output 
increases  to  410  TW,  in  excess  of  what  tungsten  wire  arrays  have  achieved  on  Z.  Figures  5-7  show 
that  carbon  has  the  capability  to  radiate  as  a  blackbody  at  200  eV  from  the  K  shell  but  only  for 
masses  at  or  in  excess  of  20  mg/cm. 

Because  the  size  of  the  plasma  remained  fixed  at  4  mm  diameter  in  all  of  the  Figure  3-7 
calculations,  the  density  of  the  plasma  increased  in  proportion  to  the  mass.  One  might  question 
whether  a  density  of  8  x  1021  ions/cm3  is  attainable  in  a  dynamic  pinch,  since  such  high  densities 
have  not  been  attained  to  date.  However,  several  factors  work  in  favor  of  this  possibility.  One, 
experimentally,  a  gapless  cylinder-like  plasma  has  never  been  imploded  and  compressed  on  axis 
before.  Two,  carbon  plasmas  do  not  need  to  be  heated  beyond  200  eV.  And  three,  the  number 
of  electrons  per  ion  is  much  smaller  in  carbon  than  in  higher  atomic  number  loads.  Thus,  if  the 
same  pressures  are  generated  in  a  carbon  plasma  on  Z  as  have  been  generated  in  higher  atomic 
number  loads,  then  higher  densities  should  be  achievable  in  carbon  than  have  been  achieved  in  other 
loads.  Most  importantly,  experiments  on  plastic  loads  should  enlighten  the  Z-pinch  community 
considerably  on  the  radiative-hydrodynamics  of  Z-pinch  loads. 


Memorandum  Report  6721-93-7347  (1993). 
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II.  ANALYZING  SPATIALLY  RESOLVED  LINE  SPECTRA 

EMITTED  BY  Z  PINCHES 


SUMMARY 

The  spatial  profiles  of  key  physical  variables  such  as  electron  temperature  and  density 
are  crucial  in  determining  the  overall  radiating  efficiency  of  a  Z  pinch  as  well  as  the  emitted 
photon  energy  distribution.  The  work  reported  in  this  chapter  uses  spatially  resolved  x-ray 
lines  emitted  from  shots  taken  on  the  4  MA  Double  EAGLE  generator  [C.  Deeney  et  al., 
Phys.  Rev.  A  44, 6762  (1991)]  at  Maxwell  Physics  International  to  infer  such  profiles.  The 
loads  were  arrays  of  15  wires  consisting  of  95%  Al,  5%  Si,  producing  nearly  50  kJ  of  K-shell 
radiation  with  relatively  long  pulsewidths  of  41-43  ns.  We  find  that  the  spatial  profile  of 
the  He-like  intercombination  to  resonance  line  ratio  is  critically  affected  by  self-absorption 
and  re-emission  (scattering)  of  resonance  line  photons  emitted  from  the  hot,  dense  interior 
of  the  pinch.  This  ratio  cannot  be  correctly  interpreted  without  taking  the  scattering  into 
account.  Also,  at  densities  somewhat  lower  than  those  achieved  by  these  pinches,  this  line 
ratio  becomes  more  temperature  than  density  dependent.  The  inferred  temperature  gradient 
is  much  less  steep  than  that  which  had  been  diagnosed  previously  for  a  pinch  emitting  a 
much  shorter  x-ray  pulsewidth  of  13  ns. 

I.  INTRODUCTION 

When  pulsed-power  generators  implode  well-designed  loads  with  currents  of  several 
MA  or  more,  total  x-ray  yields  of  100  kJ  or  more  are  typically  produced  from  the  resulting 
Z  pinch1-13.  The  highest  energy  x-rays  are  often  derived  from  emission  processes  in  the 
H-  and  He-like  ionization  stages  of  the  load  element(s),  i.e.,  the  K  shell.  X-ray  yields  from 
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the  K  shell  are  correlated  with  the  temperature  and  density  profiles  of  the  pinch.  It  is  these 
profiles  which  determine  how  much  of  the  imploded  load  is  both  heated  to  temperatures  and 
compressed  to  densities  sufficient  for  the  K-shell  ions  to  exist  in  abundance  and  to  radiate 
efficiently.  Accordingly,  major  goals  of  analyzing  experimental  data  are  to  determine  the 
fraction  of  the  original  load  mass  which  is  ionized  to  the  K  shell,  and,  if  possible,  to  make 
quantitative  inferences  of  the  temperature  and  density  gradients.  To  date,  much  of  the  work 
in  this  direction  has  employed  a  single  temperature-dependent  line  ratio,  usually  that  of 
the  H-like  to  He-like  a  lines,  together  with  the  measured  diameter  and  K-shell  power  of 
the  pinch,  to  deduce  an  average  temperature  and  density3’4’6’9,11’13-15.  Although  this 
approach  averages  over  the  spatial  variations  and  thus  only  indirectly  considers  the  gradients 
present,  it  has  led  to  some  significant  insights.  For  instance,  using  lighter  loads  generally 
increases  the  electron  temperature  of  the  assembled  pinch  while  reducing  its  density4. 
Also,  the  fraction  of  the  original  load  mass  which  participates  in  the  K-shell  radiation  (as 
obtained  from  the  spatially  averaged  analysis)  has  been  found  to  correlate  directly  with  the 
observed  K-shell  yields11.  Rarely  does  the  radiating  mass  fraction  approach  100%  of  the 
load,  indeed,  excellent  yielding  shots  commonly  exhibit  mass  participations4,9 ’11,13  of 
10-60%,  depending  on  the  load  element,  current,  and  load  design.  Such  low  mass  fractions 
suggest  that  there  is  room  for  further  improvement  in  the  K-shell  radiating  efficiency  of  Z 
pinches.  Clearly,  more  detailed  knowledge  of  how  the  gradients  present  in  Z  pinches  are 
formed  is  needed  along  with  the  factors  which  may  influence  their  formation. 

As  would  be  expected,  inferring  temperature  and  density  profiles  from  spectral  data 
requires  more  extensive  data  and  analysis  than  is  needed  to  infer  average  values.  The  use 
of  spatially  resolved  spectroscopy10,16  is  the  method  chosen  in  the  present  work.  Other 
techniques  which  have  been  employed  are:  the  replication  of  several  line  intensities  in 
spatially-integrated  spectra17,18  whose  differing  optical  depths  make  possible  the  inference 
of  plasma  conditions  at  various  depths  in  the  pinch,  laser  interferometry  ,  and 

on  no 

measuring  the  absorption  of  a  backlighting  laser  or  x-ray  source  ’  .  The  present  work 
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explores  the  potential  of  spatially  resolved  line  intensities  and  ratios  to  diagnose  critical 
features  of  the  gradient  structure  of  Z  pinches.  The  experiments  and  data  are  described  in 
Sec  II,  the  analysis  in  Sec.  Ill,  and  the  results  are  put  into  perspective  in  Sec.  IV. 

II.  EXPERIMENTS  AND  DATA 

The  experiments  analyzed  in  the  present  work  were  conducted  on  the  Double  EAGLE 
generator  at  Maxwell  Physics  International6.  This  pulsed  power  waterline  machine 
implodes  a  typical  wire  array  load  with  peak  currents  near  4  MA  reached  in  a  nsetime 
of  about  110  ns  for  a  short  circuit  load.  The  two  shots  considered  here,  4116  and  4117, 
employed  arrays  of  15  A1  wires  coated  with  Si  for  an  overall  mass  ratio  of  95%  Al,  5% 
Si,  and  mass  load  of  0.22  mg/cm.  The  length  of  the  arrays  was  2  cm,  and  their  initial 
diameter,  2.5  cm.  On  each  of  the  shots,  an  extensive  suite  of  electrical  and  radiation 
diagnostics  was  fielded,  including  filtered  x-ray  diodes  (XRDs),  photoconducting  detectors 
(PCDs),  calorimeters,  x-ray  spectrometers,  Rogowski  coils  to  measure  the  load  current, 
and  time-gated  pinhole  cameras  filtered  to  image  K-shell  radiation.  These  instruments 
produced  excellent  agreement  (better  than  10%)  in  all  of  the  major  electrical  and  radiation 
measurements  when  comparing  the  two  shots.  Therefore,  we  have  analyzed  the  specific 
results  from  shot  4116  as  representative  of  the  experiments.  A  K-shell  yield  of  47  kJ 
was  produced  within  a  full  width  at  half  maximum  (FWHM)  of  43  ns.  Despite  this  quite 
long  pulsewidth  (often  seen  when  relatively  few  wires  are  employed),  the  K-shell  yield  of 
47  kJ  is  excellent,  not  far  shy  of  the  record  of  58  kJ  for  comparable  currents  and  photon 
energies6.  This  shot  produced  a  total  x-ray  yield  of  170  kJ.  The  current  at  peak  x-ray 
power  was  4.0  MA.  While  large  numbers  of  finer  wires  are  evidently  necessary  to  obtain 
sub-10  ns  pulsewidths  and  high  radiated  powers8,11  13,  excellent  yields,  albeit  with  long 
pulsewidths  and  low  powers,  can  still  be  attained  with  a  small  number  of  thicker  wires.  The 
present  analysis  can  contribute  to  understanding  the  differences  between  short-pulse  and 
long-pulse  pinches  as  well  as  to  the  development  and  implementation  of  spatially  resolved 
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spectroscopic  diagnostics  for  Z  pinches  in  general. 

The  instrument  upon  which  much  of  the  foregoing  analysis  is  based  is  a  high-resolution, 
radially  resolving  Johann  spectrometer.  By  aligning  the  slit  along  the  pinch  axis,  radially 
resolved,  but  axially  and  time-integrated  profiles  of  individual  x-ray  lines  were  obtained. 
Typical  spatial  resolution  was  ~  0.5  mm,  and  spectral  resolution  (A/AA)  ~  1.25  x 
103.  Emitted  by  only  5%  of  the  plasma,  the  Si  lines  are  much  less  opaque  than  their  A1 
counterparts  and  have  provided  most  of  the  spectroscopic  data  analyzed  in  the  present  work. 
The  He-like  Si  Is  -ls2p  P  intercombination  (IC)  line  is  optically  thin  and  its  Doppler- 
broadened  profile  is  consistent  with  an  ion  temperature  of  5.4  keV.  Strictly  speaking,  this 
time-integrated  data  can  only  provide  an  upper  limit  for  the  ion  temperature.  However, 
given  a  nearly  stationary  pinch  during  the  43  ns  when  about  85%  of  the  K-shell  radiation  is 
emitted,  this  measurement  is  likely  not  compromised  by  bulk  motional  broadening. 

For  shot  41 16,  the  measured  spatial  profiles  of  the  He-like  Si  IC  to  He-a  resonance  line 
(ls2-ls2p1P)  ratio,  the  (He-a+IC)  intensity,  and  the  ratio  of  the  H-like  Si  ls-2p  (Ly-cr)  to 
(Hecr+IC)  line  are  shown  in  Fig.  1  (solid  lines).  Note  the  gradual  and  monotonic  decline 
with  radius  of  all  three  of  these  quantities.  The  intensity  profile  of  the  Si  He-a+IC  lines 
indicates  a  spatial  FWHM  of  3.5  mm  averaged  over  the  length  of  the  pinch.  The  radial 
profile  of  the  He-like  A1  ls2-ls3p1P  (He-/3)  line  (not  shown)  has  a  FWHM  of  3.8  mm. 
This  indicates  that  any  spatial  separation  of  the  A1  and  Si  is  mild,  and  in  the  analysis  we 
will  assume  that  the  A1  and  Si  are  fully  intermixed.  The  overall  experimental  uncertainty 
of  the  quantities  plotted  in  Fig.  1  is  estimated  to  be  15%.  Beyond  a  radius  of  2.5  mm, 
some  radiation  is  still  detected,  especially  in  He- a,  but  the  signal  is  so  close  to  background 
that  the  uncertainty  quickly  rises  to  the  30-60%  level.  Therefore,  we  have  confined  our 
quantitative  analysis  to  within  2.5  mm  of  the  pinch  axis. 

III.  ANALYSIS 
A.  Model  and  fitting  procedure 
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The  principal  analytical  tool  which  we  have  employed  is  an  ionization  dynamics  and 
radiation  transport  model  of  a  one-dimensional  (ID)  cylindrical  plasma  model  containing 
detailed  configuration  atomic  structure  for  the  Li-,  He-,  and  H-like  ionic  stages  of  A1  and 
Si.  This  model  is  described  in  more  detail  in  Ref.  15.  It  includes  all  the  important  atomic 
collisional  and  radiative  processes  coupling  the  levels.  To  get  accurate  photoexcitation  rates 
and  to  correctly  calculate  the  self-  absorption  and  power  outputs,  all  optically  thick  lines  are 
transported,  as  well  as  the  two  free-bound  continua  arising  from  radiative  recombination 
to  the  K-shell  ground  states.  The  model  solves  for  the  steady-state  populations  of  the 
individual  levels  assuming  collisional-radiative  equilibrium  (CRE).  Spatial  profiles  of  the 
spectral  lines  are  also  calculated.  Once  the  populations  are  calculated,  the  (monodirectional) 
intensity  I  of  an  x-ray  transition  emergent  along  a  chord  perpendicular  to  the  pinch  axis 
is  given  by  I  =  (Aul  /  4n)hv  fchord  Nu(s)Pe(s)ds  where  Au,  is  the  spontaneous 
decay  rate  for  the  transition,  Nu  is  the  upper  level  population  density  as  a  function  of  the 
distance  s  along  the  chord,  and  Pe  (s)  is  the  line  escape  probability  for  a  Voigt  profile  . 
This  aspect  of  the  model  calculation  is  carried  out  for  41  chords  whose  impact  parameters 
range  from  r=0  (intersecting  the  axis)  to  within  1.5%  of  the  outer  radius  of  the  pinch. 

Use  of  the  calculated  K-shell  spectrum  to  fit  the  data  (intensities,  line  ratios,  powers, 
spatial  profiles)  is  similar  to  that  described  in  Ref.  18.  The  assumed  temperature  and 
density  profiles  are  iterated,  typically  requiring  dozens  to  hundreds  of  separate  calculations 
to  find  the  best  fit.  Initially,  simple  but  versatile  analytic  forms  for  the  profiles  are  adopted, 
such  as  power  laws  or  generalized  Gaussians.  As  the  calculations  approach  the  data,  hand 
adjustments  of  the  profiles  based  on  the  observed  variable  sensitivities  are  made  to  fine 
tune  the  spectrum  and  zero  in  on  an  optimum  fit.  Based  on  the  number  of  calculations, 
and  the  fact  that  similar  but  not  identical  profiles  degrade  the  agreement  with  the  data,  we 
believe  that  the  fits  are  unique  in  the  sense  that  no  profile  differing  significantly  from  the 
one  deduced  could  fit  the  data  as  well.  However,  they  cannot  be  proven  to  be  so. 

B.  The  intercombination  to  resonance  line  ratio  diagnostic 
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The  theory  behind  the  use  of  the  He-like  IC  to  resonance  line  ratio  to  diagnose  a  plasma’s 

o  e 

average  density  has  been  reviewed  by  Griem  ,  who  also  provides  some  early  references. 
Some  recent  uses  of  this  line  ratio  to  diagnose  densities  of  plasma  focus  discharges  and 
laser-driven  plasmas  are  described  in  Refs.  26  and  27,  respectively.  The  density  sensitivity 
of  this  ratio  arises  from  the  much  larger  spontaneous  decay  rate  of  the  resonance  line  (for  Si, 
3.76  x  1013  sec-1),  compared  to  that  of  the  intercombination  line  (1.58  x  1011  sec-1).  As 
the  electron  density  increases,  the  intercombination  transition  is  increasingly  quenched  by 
collisional  de-excitation  of  the  ls2p  P  level,  whereas  the  resonance  line  continues  to  radiate 
freely  because  its  very  high  spontaneous  decay  rate  well  exceeds  that  of  electron  collisional 
de-excitation.  Therefore,  as  the  density  increases,  the  ratio  of  the  intercombination  to 
resonance  lines  decreases.  This  diagnostic  principle  is  illustrated  in  Fig.  2  where  the 
results  of  CRE  calculations  for  a  3  mm  diameter  cylinder  are  presented  in  the  form  of 
isocontours  of  this  (Si)  ratio  vs.  electron  temperature  and  Si  ion  density  for  a  95%  Al, 
5%  Si  plasma.  Above  ion  densities  of  ~  10  cm  ,  the  contours  are  mostly  vertical, 
indicating  the  dominance  of  density  dependence,  and  do  decrease  in  value  with  increasing 
density,  consistent  with  the  theory  summarized  above.  This  density  range  includes  those  of 
most  Z  pinches,  therefore,  this  line  ratio  can  still  be  considered  to  be  an  effective  density 
diagnostic  for  such  plasmas.  However,  below  this  density,  the  contours  become  horizontal, 
indicating  an  almost  pure  temperature  dependence  at  Si  ion  densities  lower  than  ~  10 
cm-3  (electron  densities  lower  than  ~  2  x  1018  cm-3).  What  accounts  for  this  unexpected 
result  showing  that  this  line  ratio  can  be  used  for  temperature  diagnosis  for  low  density 
plasmas? 

The  answer  to  this  question  lies  in  the  collisional  excitation  rates  from  the  He-like  Is2 
1S  ground  state  to  the  n=2  singlet  and  triplet  levels  from  which  arise  the  resonance  and 
intercombination  lines,  respectively.  Our  calculations  show  that,  as  function  of  temperature, 
excitation  to  the  singlet  level  rises  much  more  steeply  than  that  to  the  triplet.  For  instance, 
as  the  electron  temperature  increases  from  0.5  to  1.5  keV,  the  excitation  rate  to  the  Si  triplet 
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level  increases  by  a  factor  of  5,  but  to  the  singlet  level,  by  a  factor  of  15.  The  calculations 
of  Zhang  and  Sampson28  demonstrate  the  same  trend.  Quantum  mechanically,  it  is  the 
differing  effects  of  exchange  as  a  function  of  projectile  energy  on  the  direct  and  spinflip 
transitions  which  produce  the  relative  increase  in  excitation  of  the  resonance  line.  The 
temperature  dependence  shown  in  the  lower  density  sector  of  Fig.  2  follows  this  trend, 
i.e.,  as  the  temperature  increases,  the  resonance  line  intensity  increases  relative  to  that  of 
the  intercombination  line,  reducing  the  ratio.  In  this  low  density  regime,  neither  line  is 
collisionally  quenched.  At  higher  density,  the  quenching  effects  predominate,  producing 
the  usually  quoted  density  dependence  described  above. 

C.  Results  of  data  fit 

The  temperature  and  density  profiles  which  are  employed  to  reproduce  the  line  ratios 
and  intensities  of  the  data  of  Fig.  1  are  shown  in  Fig.  3.  Note  from  the  dotted  lines  of 
Fig.  1  that  both  the  Si  line  ratios  and  intensity  profiles  predicted  by  the  model  lie  within  the 
estimated  experimental  uncertainty  of  15%.  The  electron  density  profile,  when  translated 
into  Si  and  A1  ion  densities,  accounts  for  only  40%  of  the  initial  mass  load.  The  gentle 
decline  in  electron  temperature  shown  in  Fig.  3  (from  1.24  to  1.04  keV,  axis  to  outer 
edge)  is  required  to  replicate  the  gradual  decline  in  the  measured  temperature-sensitive 
Ly-o;/(He-a+IC)  ratio  shown  in  Fig.  1.  This  analysis  can  be  used  to  determine  only  the 
inner  K-shell  core  conditions  of  the  pinch,  as  a  cooler  outer  halo  exists  which  emits  only  the 
softer  L-shell  component  of  the  total  x-ray  output.  The  decline  in  temperature  and  density 
with  increasing  radius  is  also  consistent  with  the  observed  decrease  of  the  intensity  of  the 
(He-a+IC)  lines.  What  is  initially  puzzling  and  apparently  paradoxical  is  the  decrease  in 
density  with  radius,  despite  the  concurrent  decrease  of  the  IC/He-o:  ratio.  According  to 
the  standard  theory  outlined  above  in  Sec.  II  B  and  Fig.  2,  a  decreasing  line  ratio  should 
indicate  increasing  density. 

The  explanation  of  this  surprising  result  is  the  presence  of  line  photon  scattering 
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combined  with  the  drastically  differing  absorption  cross  sections  of  the  resonance  and 
intercombination  lines.  The  upper  levels  of  both  lines  are  competitively  populated  by 
collisional  excitation,  and  also  by  radiative  excitation  by  the  flux  of  line  photons  emitted 
from  the  hotter,  denser  pinch  interior.  The  collisional  excitation  coefficients  of  both  lines 
drop  sharply  with  decreasing  electron  temperature,  and  in  the  outer  regions  of  the  K-shell 
emitting  core,  radiative  excitation  becomes  increasingly  dominant  for  the  resonance  line 
as  is  shown  in  Fig.  4.  Note  that,  at  the  outer  edge,  the  He-o;  resonance  line  is  populated 
twice  as  often  by  photoexcitation  as  by  collisional  excitation.  However,  as  also  shown 

q 

in  Fig.  4,  radiative  excitation  of  the  intercombination  line’s  P  upper  level  is  at  most  a 
1%  effect.  This  is  a  straightforward  consequence  of  the  fact  that  the  photoexcitation  cross 
section  of  the  intercombination  line  is  less  than  that  of  the  resonance  line  by  a  factor  of  230. 
With  increasing  radius,  resonance  line  emission  is  increasingly  photo-excited  by  the  flux  of 
x-rays  from  the  hotter,  denser  interior  whereas  this  population  mechanism  is  not  available 
to  the  intercombination  line.  Thus,  only  the  spatially  resolved  intensity  of  the  resonance 
line  is  enhanced  by  this  scattering  mechanism,  reducing  the  IC/resonance  ratio  even  though 
the  electron  density  is  decreasing. 

Failure  to  properly  take  this  mechanism  into  account  can  lead  to  erroneous  and  seriously 
misleading  density  inferences.  Fig.  5  shows  the  electron  density  profile  that  is  derived  from 
the  observed  IC/resonance  line  ratio  when  x-ray  scattering  is  neglected.  Note  that  the 
density  increases  with  radius,  in  accord  with  the  usual  interpretation  of  this  line  ratio. 
When  this  new  density  profile  is  used,  the  spatially  resolved  intensity  of  the  (Hect+IC)  lines 
is  found  to  increase  with  increasing  radius,  as  also  plotted  in  Fig.  5,  and  in  contrast  to  the 
data. 

The  reversal  of  the  normal  density  dependence  of  the  IC/resonance  line  ratio  is  a 
space  dependent  phenomenon.  However,  most  of  the  energy  for  driving  both  lines  comes 
from  collisional  excitation  in  the  hot,  dense  interior  near  the  axis.  Therefore,  the  spatially 
integrated  ratio  would  preserve  the  usual  density  dependence.  In  Fig.  6,  the  calculated 


18 


spatial  profiles  of  the  IC/resonance  line  ratio  are  displayed  for  three  different  density  profiles, 
the  inferred  profile  of  Fig.  3  for  Double  EAGLE  shot  4116,  and  also  profiles  representing  a 
doubling  and  halving  of  this  original  density  distribution  at  each  point.  Note  that  for  each 
profile,  the  ratio  declines  with  radius  along  with  the  electron  density.  However,  for  the 
highest  electron  density  the  overall  ratio  is  lowest  and  vice-versa,  in  accord  with  the  usual 
interpretation  of  its  density  dependence.  Therefore,  in  a  spatially  unresolved  spectrum,  a 
high  IC/resonance  ratio  does  indicate  a  relatively  low  average  density  and  a  low  ratio  a  high 
average  density.  But,  to  correctly  infer  density  profiles  within  this  overall  density  range,  it 
is  necessary  to  account  for  photoexcitation  and  its  spatial  dependence. 

How  do  the  averages  of  spatial  temperature  and  density  profiles  compare  with  those 
that  would  be  derived  from  the  more  often  employed  spatially  integrated  spectroscopic 
analysis  using  the  dominant  pinch  component,  i.e.,  Al?  As  described  in  Sec.  I,  the  spatially 
averaged  analysis  uses  the  measured  K-shell  power  and  pinch  diameter  together  with  a 
temperature-dependent  line  ratio  such  as  Ly-o;/(He-Q:+IC)  to  infer  an  electron  temperature 
and  ion  density.  For  this  shot,  the  time-integrated  Ly-a/(He-a+IC)  ratio  is  1.3,  the  peak 
Al  K-shell  power  is  0.5  TW/cm,  and  the  diameter  of  the  pinch  in  the  Al  He-/3  line  is  3.8 
mm.  Employing  CRE-calculated  contour  plots  similar  to  those  presented  in  Ref.  15,  and 
using  the.  peak  K-shell  power  quoted  above,  an  Al  ion  density  of  1.4  x  10  cm  and 
electron  temperature  of  0.98  keV  are  obtained.  These  values  are  in  only  fair  agreement 
with  the  density  of  8.0  x  1018  cm-3  and  temperature  of  1.15  keV  which  are  obtained 
by  spatially  averaging  (over  3.8  mm)  the  profiles  of  Fig.  3.  However,  when  the  0.25 
TW/cm  average  K-shell  power  during  the  pulse,  (excluding  the  100  ns  “foot”  which  begins 
when  the  K-shell  power  has  declined  to  20%  of  its  peak  value)  is  employed  as  input  to  the 
Al  spectroscopic  analysis  instead  of  the  peak  power,  much  better  agreement  is  achieved, 
namely  an  Al  ion  density  of  9.8  x  1018,  and  electron  temperature  of  1.11  keV.  Therefore 
the  two  methods  are  complementary  and  consistent  provided  that  the  time-averaged  power 
is  used  for  the  spatially  integrated  analysis  to  correspond  to  time  integrated  intensities  and 


19 


line  ratios  employed  in  the  spatially  resolved  analysis. 

IV.  CONCLUDING  REMARKS 

Much  of  the  spectroscopic  analysis  of  Z-pinches  that  has  been  carried  out  to  date 
has  employed  spatially  integrated  data.  Such  a  comparatively  simple  approach  implicitly 
averages  over  the  gradients  which  are  present  to  varying  degree  in  all  pinches.  However,  it 
has  nonetheless  yielded  quantitative  connections  between  the  radiative  output  of  Z  pinches 
and  the  atomic  number  of  the  load  and  the  currents  and  pulsewidths  with  which  they  are 
driven.  Thus,  even  for  shots  in  which  the  best  K-shell  yields  are  obtained,  it  is  common 
for  well  under  half  the  load  mass  to  be  producing  the  K-shell  radiation.  This  presents  the 
opportunity  to  enhance  the  K-shell  yield  by  designing  experiments  that  involve  more  of  the 
mass  heated  to  K-shell  temperatures.  However,  to  effectively  carry  out  such  experiments, 
the  gradients  which  hold  the  key  to  the  lack  of  full  mass  partication  must  be  routinely 
diagnosed. 

The  work  described  in  this  report  shows  one  example  of  how  this  diagnostic  challenge 
can  be  met  through  the  use  of  detailed  spatially  resolved  spectra.  Such  data  were  analyzed 
from  two  identical  shots  on  the  4  MA  Double  EAGLE  generator  which  produced  consistent 
K-shell  yields,  approaching  50  kJ,  from  a  mixture  of  95%  Al,  5%  Si.  In  the  analysis,  we 
noted  that  the  ratio  of  the  intercombination  to  He-a  lines,  normally  considered  to  be  a 
density  diagnostic,  could  be  used  as  an  excellent  temperature  diagnostic  for  densities  up  to 
two  orders  of  magnitude  lower  than  those  normally  encountered  in  Z  pinches.  For  Z-pinch 
densities,  this  ratio  can  be  an  effective  density  diagnostic,  and  we  have  reconciled  its  use 
with  that  of  absolute  powers  to  determine  the  average  and  spatially  resolved  densities.  In 
diagnosing  the  density  profile  within  the  pinch  from  spatially  resolved  data,  however,  it 
is  necessary  to  properly  take  into  account  photoexcitation  to  avoid  seriously  misleading 
inferences  (see  Fig.  5). 

To  fully  realize  the  benefits  of  this  kind  of  data  analysis  a  greater  variety  of  shots  needs 
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to  be  systematically  analyzed  for  their  gradient  structure.  However,  it  is  interesting  to 
compare  the  present  results  to  those  previously  obtained  from  shot  70  on  the  Z  generator 
at  Sandia  National  Laboratories18,  which  was  also  predominantly  Al.  From  the  Z  shot,  a 
K-shell  yield  of  122  kJ  was  obtained  within  a  13  ns  pulse,  much  shorter  than  the  41-43  ns 
characteristic  of  the  presently  considered  Double  EAGLE  experiments.  The  K-shell  yield 
on  Z  was  far  below  the  implosion  kinetic  energy  of  ~  1.1  MJ,  whereas  on  Double  EAGLE, 
the  K-shell  yield  of  47  kJ  was  much  closer  to  the  kinetic  energy  of  ~  63  kJ,  with  the  total 
radiative  yield  of  170  kJ  nearly  triple  the  coupled  kinetic  energy.  On  Z,  the  total  yield 
and  coupled  kinetic  energy  were  about  equal  for  the  Al  shots.  The  temperature  gradient 
diagnosed  on  Z  shows  the  electron  temperature  falling  from  1 .6  keV  to  0.2  keV  within  1  mm, 
whereas  a  much  gentler  fall  in  temperature  from  1.2  to  1.0  keV  within  2.5  mm  is  inferred 
for  the  Double  EAGLE  shot.  These  sharply  contrasting  factors  in  yield  vs.  coupled  kinetic 
enegy,  temperature  gradient,  and  K-shell  pulsewidth,  may  be  related  to,  or  a  consequence 
of,  differences  in  plasma  opacity,  implosion  dynamics  or  energy  coupling  mechanisms, 
such  as  the  “magnetic  bubble”  energy  deposition  mechanism  recently  proposed  .  Further 
investigation  and  more  comprehensive  and  systematic  analysis  is  needed  to  verify  these 

relationships. 
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FIGURE  CAPTIONS 


FIG.  1.  The  indicated  experimental  and  theoretical  line  ratios  and  intensities  are  plotted 
vs.  radius  for  Double  EAGLE  shot  4116.  The  experimentally  measured  quantities  are 
indicated  by  the  triangles  connected  by  solid  lines,  the  model  fit  results  given  by  the  dashed 
line  curves.  Overall  experimental  uncertainty  is  estimated  at  15%. 

FIG.  2.  Calculated  isocontours  of  the  intercombination  to  resonance  line  ratio  for  the  He- 
like  Si  ion  are  shown  as  a  function  of  electron  temperature  and  Si  ion  density  for  a  3  mm 
diameter  cylindrical  plasma  composed  of  95%  Al,  5%  Si. 

FIG.  3.  Electron  density  and  temperature  profiles  inferred  from  the  data  presented  in  Fig. 
1,  for  Double  EAGLE  shot  41 16. 

FIG.  4.  The  ratio  of  radiative  to  collisional  excitation  is  plotted  for  the  profiles  of  Fig.  3, 
for  the  Ly-Ct,  He-a,  and  intercombination  lines  of  K-shell  Si  ions. 

FIG.  5.  The  electron  density  inferred  when  photoexcitation  is  neglected,  and  the  (He-a+IC) 
intensity  profile  predicted  from  this  inferred  density  are  plotted  vs.  displacement  from  the 
axis.  Also  shown  are  the  electron  density  obtained  from  the  full  model  as  well  as  the 
measured  (He-a+IC)  intensity  profile. 

FIG.  6.  The  calculated  spatial  profiles  of  the  IC/resonance  line  ratio  are  plotted  for  the 
electron  density  profile  of  Fig.  2  as  well  as  for  profiles  in  which  this  density  distribution 
has  been  halved  and  doubled  at  each  point  (central  electron  density  is  indicated). 
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III.  AN  EFFICIENT  NON-LTE  RADIATION  TRANSPORT  MODEL  SUITABL 
FOR  MULTIDIMENSIONAL  HYDRODYNAMICS  CALCULATIONS 


I.  Introduction 

A  straightforward  approach  to  establishing  a  theoretical  relationship  between  the  initial 
asymmetries  or  the  Rayleigh  Taylor  (RT)  instabilities  of  a  z-pinch  and  its  ability  to  efficiently 
radiate  K-shell  photons  would  be  to  perform  two-dimensional  (2D)  or  three-dimensional  (3D),  time 
dependent,  non-local  thermodynamic  equilibrium  (non-LTE),  full  radiation  transport,  magneto¬ 
hydrodynamic  (MHD)  calculations.  However,  this  capability  has  never  been  demonstrated  by  the 
z-pinch  community  because  the  computer  requirements  for  performing  the  zone-to-zone  coupling 
needed  to  solve  the  equation  of  radiative  transfer  in  a  2D  or  a  3D  plasma  are  severe.  In  the 
past,  our  2D  MHD  model,  PRISM,1  which  stands  for  Plasma  Radiating  Imploding  Source  Model, 
inadequately  modeled  the  radiation  field  by  assuming  either  that  the  plasma  was  optically  thin 
or  that  a  local  escape  model  could  be  used  to  model  some  of  the  photon  self-absorption  physics. 
Another  model  widely  used  by  the  z-pinch  community,  the  LTE  radiative  diffusion  model,  is  also 
inadequate  for  modeling  the  radiation  physics  of  plasma  radiating  source  (PRS)  loads  designed  for 
K-shell  emission  because  neither  the  diffusion  approximation  nor  the  LTE  assumption  are  valid  for 
K-shell  plasma  conditions  produced  by  present  day  Z-pinches. 

The  transport  of  radiation  within  an  optically  thick  plasma  is  both  a  nonlinear  and  a  nonlocal 
problem.  The  source  functions  and  absorption  coefficients  throughout  the  plasma  determine  the 
amount  of  radiation  that  escapes  the  plasma  as  well  as  the  local  radiation  energy  densities.  The 
source  functions,  in  turn,  depend  on  these  densities  through  the  photoexcitation  and  ionization 
processes  they  stimulate,  which  influences  the  ion  and  excited  state  population  densities  in  the 
plasma.  The  transport  problem  is  solved  exactly  only  when  the  radiation  field  and  populations  are 
determined  self-consistently.  For  these  reasons  the  optically  thin  and  radiative  diffusion  models 
are  not  viable  candidates  for  modeling  radiation  transport  in  multi-dimensional,  K-shell  emitting 
plasmas.  We  have  therefore  developed  a  model  that  retains  the  nonlocal  and  nonlinear  nature 
of  the  transport  problem  but  only  approximates  the  exact  solution  by  sacrificing  some  of  the 
self-consistency  inherent  in  an  exact  solution  to  the  transport  equation.  This  model  is  called  the 
Tabulated  Collisional-Radiative-  Equilibrium  (TCRE)  model. 
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This  model  is  described  in  Section  II  along  with  three  other  models  with  which  it  is  compared 
later  in  the  paper.  In  Section  II. a,  a  Global  Collisional-Radiative-Ecjui librium  (GCRE)  model  is 
discussed.  This  model  gives  accurate  solutions  to  the  exact  transport  problem,  and  it  therefore 
serves  as  a  benchmark  for  testing  the  TCRE  model.  Unfortunately,  the  self-consistency  demanded 
by  the  GCRE  model  makes  it  computationally  prohibitive  for  use  in  multi-dimensional  calculations 
on  most  computer  platforms  at  this  time.  In  Section  H.b,  a  less  accurate  on-the-spot  model  that 
is  self-consistent,  but  which  has  a  simplifying  treatment  of  the  nonlocality  problem  that  is  also 
employed  in  the  TCRE  model,  is  described.  The  self-consistency  demanded  by  this  model  also 
makes  it  computationally  prohibitive.  The  TCRE  model  is  described  in  Section  n.c.  This  section 
is  concluded  with  a  brief  mention  of  an  optically  thin  model  in  Section  Hd.  Results  obtained  from 
ID,  cylindrical  geometry,  z-pinch  calculations  using  the  TCRE  model  are  compared  in  Section  HI 
with  ID  radiation  transport  calculations2  incorporating  the  three  other  models  discussed  in  Section 
n.  The  TCRE  model  is  found  to  produce  results  in  reasonable  agreement  with  the  GCRE  model.  In 
Section  IV,  the  TCRE  model,  as  implemented  in  the  2D  PRISM  code,  is  employed  to  preliminarily 
investigate  argon  Z-pinch  experiments  similar  to  those  performed  on  the  Z-machine  at  Sandia 
National  Laboratories.3  These  2D  calculations  employ  the  measured  initial  gas  distribution  at  the 
cathode  (the  axial  dependence  of  the  gas  distribution  is  ignored)  and  the  K-shell  output  is  computed 
as  a  function  of  load  mass.  Comparisons  are  made  with  ID  model  results  of  these  same  experiments 
to  examine  the  role  that  the  Rayleigh  Taylor  instability  plays  in  affecting  the  production  of  K-shell 
emission.  This  work  is  summarized  in  Section  V. 

II.  Models  for  Radiation  Transport  and  kinetics 

The  radiation  field  in  a  plasma  alters  the  atomic  level  populations  through  photo-pumping 
while  the  populations  determine  the  field  through  emission  and  absorption  processes.  Therefore 
radiation  transport  and  population  kinetics  are  coupled  and  can  only  be  solved  in  a  unified  approach. 
The  transport  equation  for  the  radiation  intensity  Iv{ r0,  t,  Cl)  at  frequency  v,  at  position  r0,  at  time 
f,  traveling  in  the  unit  direction  Cl,  is 

Cl  •  V/j/  +  KVIV  =  ju  ,  (1) 

where  kv  is  the  frequency  dependent  absorption  coefficient  in  units  of  inverse  length  and  ju  is 
the  emission  coefficient  at  frequency  v  in  units  of  energy  per  unit  volume,  per  unit  time,  per  unit 


33 


frequency,  per  steridian.  The  finite  light  transit  time  across  the  medium  has  been  neglected  in 
eqn.(l),  as  is  appropriate  for  Z-pinches.  Integrating  /„  over  all  directions  at  a  canonical  position  r0 
gives  the  mean  intensity  J„(r0)  =  (1/4tt)  f  RdQ.  For  the  population  kinetics  consider  a  simple 
atomic  level  model  with  ground  level  ionization  stages  j  and  j  +  1,  and  stage  j  having  a  single 
excited  state  k.  The  number  density,  nk(r0,  f),  in  the  excited  state  k  of  the  jr’th  ionization  stage  at 
position  rc  and  time  t,  is  calculated  from 


oTljCj'k  TleTlkDkjj  "b  '^e^,kC'k,j+l  T  ’^je'^'j-{-lRj+l,k 


-nk  J  a*'j+lj^dv +  nj  J  aJu'k^dv  -  nkAkJ  . 


The  electron  collisional  excitation  rate  Cj)k  is  from  the  ground  state  of  ion  j  to  the  excited  state  k, 
and  Dk,j  is  the  corresponding  de-excitation  rate.  Three-body  recombination  from  the  next  higher 
ionization  state  is  given  by  Bj+i>k,  and  Ckj+ 1  is  the  electron  ionization  rate.  The  remaining 
four  terms  represent  radiative  processes.  Rj+i,k  is  the  radiative  recombination  rate.  Dielectric 
recombination  (not  explicitly  displayed)  can  be  included  with  the  same  functional  form  as  the  R 
term.  The  frequency  integrals  are  over  the  photo-ionization  cross-section  crk’i+l  and  the  photo¬ 
excitation  cross-section  cr£fc.  These  cross-sections  include  stimulated  emission  processes.  The 
transition  probability  Akj  accounts  for  radiative  decay  of  the  excited  state.  The  populations  from 
eqn.(2)  are  coupled  to  the  radiation  field  of  eqn.(l)  through  the  expressions  for  the  emission  and 
absorption  coefficients.  For  the  line  emitted  by  the  excited  state  k  and  absorption  by  the  ground 
state  j  one  has 

jv  =  ^nkhvkijAkJ4>k>:>  and  ku  =  nj<Tl’k  ,  (3) 

where  h  is  Planck’s  constant  and  4>k^  is  the  normalized  line  profile  function:  f  (p^df  =  1.  The 
radiative  power  density,  the  net  radiation  energy  lost  or  gained  per  unit  volume  per  unit  time  from 
each  point  in  the  plasma,  is 

A  =  47T  f  (j„  -  nuJv)du  .  (4) 


This  expression  can  be  positive  or  negative  depending  on  the  position  in  the  plasma  and  enters  as 
a  loss  or  gain  term  in  the  electron  energy  equation,  respectively.  The  integral  of  A  over  the  whole 
plasma  volume  is  the  emitted  radiated  power  from  the  plasma. 

a.  Global  Collisional  Radiative  Equilibrium 
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Four  models  are  employed  in  this  paper  to  solve  the  above  coupled  system  of  equations 
describing  the  radiation  field  and  the  atomic  level  populations.  The  assumptions  of  each  approach 
will  be  described  to  clarify  the  essential  differences.  The  first  approach  is  the  Global  Collisional 
Radiative  Equilibrium  (GCRE)  model  in  which  the  population  kinetics  of  eqn.(2)  are  solved  in  the 
steady-state.  This  model  assumes  that  the  dynamic  timescale  for  plasma  motion  is  long  compared 
to  the  timescale  for  collisional  or  radiative  processes,  which  is  accurate  for  most  K-shell  x-ray 
sources.  This  assumption  is  not  however  intrinsic  to  the  model  and  other  work  has  included  time- 
dependent  ionization  for  studies  involving  the  production  of  high  energy  x-rays.4’5  Other  than  the 
condition  of  steady-state,  the  complete  non-linear  coupling  between  the  transport  and  the  kinetics, 
which  enters  through  the  J„  terms  in  eqn.(2),  is  included  and  comprises  the  CRE  component  of 
the  model.  From  the  solution  to  eqn.(l)  the  photons  which  lead  to  photo-pumping  (ionization 
or  excitation)  can  arise  from  different  parts  of  the  plasma  than  the  position  rc.  This  leads  to 
a  non-local  spatial  coupling  between  any  two  positions,  i.e.,  zones  in  a  numerical  simulation. 
Photons  emitted  from  position  rG  and  absorbed  at  a  different  position  ? i  alter  the  populations  and 
emissivity  at  the  latter  point,  which  in  turn  affects  the  populations  at  the  originating  position  r0  due 
to  absorption  of  photons  from  rx.  This  feedback  can  be  accounted  for  by  calculating  a  matrix  of 
radiative  coupling  coefficients  connecting  each  zone  to  every  other  zone  for  each  frequency.  In  a 
full  dynamic  simulation,  the  solution  for  the  populations  and  radiation  field  at  a  particular  timestep 
is  then  iterated  while  conserving  electron  thermal  plus  excitation/ionization  energy  in  each  zone 
separately.  This  procedure  provides  a  self-consistent  solution  for  both  the  radiation  field  and  the 
level  populations. 

In  some  cases  the  matrix  of  coupling  coefficients  at  each  frequency  can  be  replaced  by  a 
similar  matrix  for  each  emission  line.  In  essence  the  set  of  frequencies  covering  a  line  is  integrated 
over  and  the  matrix  components  provide  the  probability  that  a  line  photon  emitted  from  ri  will 
transit  the  plasma  and  be  absorbed  at  r0.  Apruzese6-8  provides  frequency  and  angle  averaged 
formulas  for  the  line-profile-averaged  transit  probabilities  using  Lorentz,  Doppler,  or  Voigt  line 
profiles.  A  similar  procedure  can  be  used  for  bound-free  edges  arising  from  photo-ionization.9 
The  advantage  of  the  transit  probability  approach  over  that  of  multi-frequency  is  a  reduction  in 
the  size  of  the  matrix  of  coupling  coefficients  since  typically  10-30  frequencies  must  be  used  to 
cover  each  line.  In  either  case,  however,  the  adjective  “Global”  as  part  of  the  model  designation 
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emphasizes  the  non-local  character  of  the  solution  scheme  coupling  the  the  radiation  transport  of 
eqn.(l),  the  level  kinetics  of  eqn.(2),  and  the  emissivity  and  absorption  coefficients  of  eqn.(3).  The 
GCRE  model  has  been  used  to  calculate  the  radiative  emission  from  a  Z-pinch  implosion  in  much 
of  the  work  by  the  Radiation  Hydrodynamics  Branch  at  the  Naval  Research  Laboratory  during  the 
last  twenty  years,  see,  e.g.  references  4-10.  In  their  papers  the  approach  is  referred  to  as  CRE.  All 
of  the  dynamic  calculations  of  Z-pinch  implosions  which  employ  GCRE  have  been  limited  to  1-D 
but  include  many  excited  states  per  ionic  stage.  The  above  GCRE  approach  in  2-D,  which  requires 
a  matrix  of  radiative  coupling  coefficients  between  all  zones  at  each  emission  feature  is  at  present 
computationally  prohibitive  for  most  platforms. 

b.  “On-the-Spot” 

A  second  model  employs  the  “on-the-spot”  approximation  of  Holstein11  and  Osterbrock12.  In 
this  approximation  if  a  photon  is  absorbed  anywhere  in  the  plasma  it  is  assumed  to  be  absorbed  in 
the  emitting  zone.  The  physical  basis  of  this  approximation  is  that  most  photons  that  do  not  escape 
the  plasma  are  those  emitted  in  the  highly  opaque  cores  of  spectral  lines  and  are  therefore  absorbed 
nearby.  It  is  most  accurate  in  the  extreme  limits  of  either  a  highly  opaque  or  optically  thin  plasma. 
In  the  present  implementation  of  this  approximation  the  steady-state  condition  is  again  assumed, 
photo-ionization  in  eqn.(2)  is  neglected,  and  for  each  emission  line  one  determines  the  probability 
that  a  photon  emitted  from  position  rc  will  escape  all  the  way  to  the  edge  of  the  plasma,  i.e.,  the 
probability-of-escape  V.  Holstein11  first  developed  an  expression  for  the  mono-directional  escape 
probability  for  a  completely  frequency  redistributed  Doppler  line  profile  at  large  optical  depth.  In 
a  multi-dimensional  plasma  simulation  one  requires  an  average  of  this  escape  probability  over  all 
possible  directions.  For  the  line  emission  from  level  k  down  to  j: 

VkjVc)  =  Ay  j  tofoda .  (5) 

where  the  the  optical  depth  for  the  k-to-j  transition  r*,J-  covers  the  path  from  rc  to  the 
edge  of  the  plasma  along  the  direction  Cl.  A  ray  trace  method  is  employed  to  calculate 

=  Y2m  «£’J(rm)5sm,  where  Ssm  is  the  path  length  through  zone  m  in  the  Cl  direction  and 
KUk  .  (rm)  is  the  absorption  coefficient  in  that  zone.  The  angular  integral  in  eqn.(5)  is  calculated  in 
3-D  by  choosing  a  discrete  set  of  ordinates  with  appropriate  weights.  The  probability-of-escape  is 
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used  to  combine  the  last  two  terms  of  eqn.(2)  leading  to 

=  0  =nenjCjtk  ~  nenkDktj  +  nenerij+iBj+itk  -  nenkCk,j+ 1  /g\ 

CJb  '  ' 

+nerij+iRj+i,k  -  nkAktjVk,j  ■ 

When  the  plasma  is  optically  thick  in  a  particular  line  (V  <  1),  a  fraction  (1  -  V)  of  the  photons 
emitted  from  position  r0  are  reabsorbed  and,  consequently,  the  radiative  decay  rate  is  effectively 
reduced  to  a  fraction  Vkj  of  its  intrinsic  atomic  value.  The  computational  advantage  of  the  “on-the- 
spot”  approximation  is  the  absence  of  a  matrix  of  radiative  coupling  coefficients  since  one  needs 
only  to  compute  and  store  the  escape  probabilities  to  the  plasma  edge  and  not  to  all  zones  within  the 
plasma.  This  is  a  simplification  for  radiation  transport  in  multi-dimensional  simulations.  However, 
an  iterative  procedure  between  the  radiation  field  of  eqn.(l)  and  the  population  kinetics  in  eqn.(6) 
is  still  employed  in  this  second  model:  the  escape  probability  of  one  line  at  some  position  depends 
on  the  level  populations  and  escape  probabilities  of  other  lines  at  all  other  positions  in  the  plasma 
through  the  absorption  coefficient  k  of  eqn.(3).  Since  a  photon  either  is  absorbed  at  its  point  of 
origin  or  escapes  the  plasma,  the  expression  for  the  volumetric  radiative  loss  term  in  eqn.(4)  is 
replaced  by 

A  =  4irjkjPk,j  ) 

where  the  sum  includes  all  emission  lines. 

c.  Tabular  Collisional  Radiative 

The  third  model  retains  the  non-local  character  of  radiation  transport  inherent  in  the  first  two 
approaches,  includes  the  non-linear  coupling  between  the  radiation  field  and  the  level  populations, 
but  replaces  the  exact  self-consistency  of  a  dynamically  iterative  solution  performed  at  each  timestep 
with  an  approximate  one  needing  no  iteration.  The  main  assumption  of  this  model  is  that  the  line 
emissivity  from  a  particular  position  x0  is  equal  to  that  from  a  uniform  plasma  whose  properties 
are  given  by  the  density,  temperature,  etc.,  at  that  position.  The  size  of  this  uniform  plasma  and 
the  probability-of-escape  are  mutually  dependent.  This  assumption  is  valid  for  optically  thick  lines 
when  the  gradient  length  scale  of  the  plasma  parameters  is  large  compared  to  the  absorption  length 
of  a  photon.  For  thin  lines  the  gradient  is  unimportant  and  the  photon  escapes  along  any  line  of 
sight.  In  this  model  the  transport  and  kinetics  are  solved  using  a  set  of  pre-calculated  tables  from 
which  jk,j  and  Kjik  are  obtained. 
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To  illustrate  how  the  tables  are  established,  consider  a  uniform  1-D  plasma  slab  specified  by 
its  thickness,  5x,  total  ion  density,  rii  =  rij  +  rik,  and  the  internal  energy,  E.  Equations  (1), 
(2),  and  (3)  are  solved  self-consistently  for  a  single  zone  in  this  ideal  plasma  geometry  by  adjusting 
the  electron  temperature  Te  and  level  populations  until  E  =  (3/2 )nekBTe  +  Yjnjej  +  J2nkek, 
where  ne  is  the  electron  density  determined  from  the  ion  densities  rij,  and  £j  (e^)  is  the  ionization 
(excitation)  energy  per  particle.  For  each  transported  line  jkj,  there  is  a  unique  relation  between 
5x  and  the  probability-of-escape  Vk,j  from  the  slab.  This  facilitates  the  construction  of  Radiative 
Equation  of  State  (REOS)  tables  of  electron  temperature  Te,  mean  charge  state  Z,  emissivity  jkj, 
and  absorption  coefficient  k jtk,  as  a  function  of  ( rii,E,Vk,j )•  In  practice  at  each  timestep  of 
a  simulation  the  magnetohydrodynamic  component  of  a  code  supplies  n*  and  E  for  every  zone. 
Then  the  probability-of-escape  is  computed  from  eqn.(5)  and  the  REOS  Tables  are  searched  for 
the  corresponding  values  of  (Te,  Z,jk,j ,  Kj,k)  for  each  transported  line  from  every  zone  whenever 
the  radiation  field  is  updated.  In  this  manner  the  tabulated  Vk,j  are  matched  to  the  probabilities-of- 
escape  for  the  entire  plasma,  not  the  local  numerical  zone.  The  procedure  is  summarized  in  Fig.l. 
We  designate  this  approach  as  the  Tabular  Collisional  Radiative  Equilibrium  (TORE)  model. 

The  inconsistency  which  arises  in  the  TORE  model  stems  from  the  use  of  tabular  data  based 
on  a  uniform  medium  in  place  of  dynamically  coupling  the  transport  and  kinetics  by  an  iterative 
approach  as  described  in  the  GCRE  and  “on-the-spot”  models  above.  This  inconsistency  also  leads 
to  a  temporal  offset  in  the  spatial  coupling  of  populations  between  distant  zones  resulting  from 
the  time  interval  between  radiation  calls  calculating  the  escape  probabilities.  The  radiation  field  is 
typically  updated  on  a  much  longer  time  scale  than  the  hydrodynamic  time  scale.  An  inconsistency 
due  to  the  non-iterative  tabular  look-up  might  be  expected  to  arise  from  two  different  transported 
lines  in  different  ionization  stages,  which  have  different  probabilities-of-escape.  This  approach 
might  give  different  values  of  Te  and  Z  in  the  same  zone  based  on  the  REOS  table  look-up  . 
However,  it  is  found  that  in  the  absence  of  bound-free  opacity,  Te  and  Z  are  for  all  practical 
purposes  dependent  primarily  on  rii  and  E.  The  weak  dependence  of  Te  and  2  on  Vk,j  that  does 
exist  is  largest  at  low  ion  density  and  negligible  at  higher  ion  density.  This  point  is  illustrated  in 
Figs.2  and  3,  where  Te  and  Z  are  plotted  as  functions  of  E  and  5x  for  an  argon  plasma  with  ion 
densities  of  7.0  x  1017  and  7.6  x  1020  cm-3,  respectively.  The  consequences  of  not  iterating  are 
described  in  the  results  section. 
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In  addition  to  the  uniform  plasma  and  “on-the-spot”  approximations  of  theTCRE  model  there 
are  other  features  worth  noting.  First,  only  one  excited  level  in  each  ionization  stage  is  included  in 
the  1-D  REOS  tables  as  a  representative  line  and  is  transported  in  the  multi-dimensional  plasma. 
This  level  is  chosen  to  correspond  to  the  strongest,  most  opaque  resonance  line.  The  emissivities 
of  other  lines  emanating  from  that  ionization  stage  are  scaled  from  the  resonant  line  using  the 
uniform  plasma  approximation  and  are  added  to  the  emissivity  of  the  representative  line.  Secondly, 
at  present  only  line  opacity  is  considered  in  the  escape  probabilities.  Finally,  opacity  effects  are 
included  in  the  bremsstrahlung  radiation  emissivity  but  not  for  the  recombination  (free-bound) 

emissivity. 

In  summary,  the  TCRE  model  retains  the  advantages  of  the  escape  probability  method 
developed  for  1-D  plasmas,  and,  as  we  will  show,  provides  a  reasonably  accurate,  computationally 
efficient,  multi-dimensional  method  for  coupling  radiation  transport  and  atomic  level  kinetics. 

d.  Optically  Thin 

A  fourth  model  is  considered  in  the  paper  for  comparison  purposes.  It  treats  the  radiation  as 
optically  thin  and  ignores  transport  effects  altogether.  This  model  simply  but  incorrectly  assumes 
that  all  created  photons  escape  the  plasma  and,  hence,  the  V  term  in  eqn.(6)  is  set  to  unity.  It  is 
included  to  demonstrate  the  extreme  differences  which  arise  in  dynamical  simulations  of  Z-pinches 
employing  this  assumption  as  compared  to  any  other  treatment  of  transport. 

III.  Model  Comparisons 

The  accuracy  of  the  TCRE  model  is  assessed  by  comparing  the  results  of  ID  Lagrangian 
MHD  calculations  that  employ  this  model  against  results  obtained  using  the  GCRE  model.  Both 
models  use  the  same  ID,  cylindrical  geometry,  ray  trace  scheme.7  The  K-shell  radiative  yield,  the 
total  radiative  yield,  stagnation  profiles  for  m  and  Tc,  and  computation  times  are  compared  for 
several  typical  Z-pinch  calculations.  These  include  150  pg/ cm  and  500  pg/cm,  uniform  fill,  3.5 
cm  radius,  argon  loads  that  are  imploded  using  a  linearly  ramped  current  profile  that  is  terminated 
when  the  energy-per-ion  coupled  to  the  load  exceeds  an  77*  value  of  2.  77*  is  the  coupled  energy- 
per-ion  divided  by  the  minimum  energy-per-ion  needed  to  ionize  to  the  K  shell.13  For  argon  this 
minimum  energy  is  40  keV  per-ion.14  In  addition,  comparisons  are  made  for  1125  pg/cm  and 
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1857  pg/cm,  2  cm  long,  titanium  wire  array  loads  that  have  initial  array  diameters  of  4.5  cm  and 
3.5  cm,  respectively.  These  loads  are  imploded  using  a  circuit  model  for  the  Z  machine  at  Sandia 
National  Laboratories  and  have  rj*  values  around  2. 

Table  I  compares  the  total  and  K-shell  yields  as  calculated  for  the  above  loads  using  four 
different  radiation  models,  the  models  are:  (1)  TCRE,  (2)  “ on-the-spot ”,  (3)  GCRE,  and  (4)  the 
optically  thin.  The  ID  calculations  are  identical  except  for  the  radiation  model.  In  comparison  to 
the  other  models,  the  optically  thin  model  predicts  too  much  radiative  energy  escaping  the  z-pinch. 
For  the  titanium  loads,  the  masses  are  large  enough  that  the  optically  thin  model  radiates  the  energy 
so  efficiently  that  the  resulting  overcooled  plasma  can  barely  ionize  into  and  radiate  from  the  K 
shell.  This  comparison  emphasizes  the  importance  of  photon  trapping. 

Since  the  TCRE  model  is  inherently  the  “ on-the-spot ’  model  with  the  exception  that  the 
uniform  plasma  approximation  is  also  implemented,  a  comparison  between  the  TCRE  and  “ on- 
the-spot' ’  models  is  a  good  indicator  of  how  much  the  uniform  plasma  approximation  is  adversely 
affecting  the  TCRE  calculation.  In  terms  of  K-shell  and  total  yields,  Table  I  shows  that  for  these 
typical  z-pinch  loads  the  uniform  plasma  approximation  is  very  good. 

The  Table  I  comparison  of  the  TCRE,  “ on-the-spot ’,  and  GCRE  models  shows  that,  with  the 
exception  of  the  K-shell  emission  of  the  1857  pg/cm  titanium  load,  the  models  are  in  good  to 
excellent  agreement.  One  reason  why  model  deficiencies  are  likely  to  show  up  for  a  load  such  as  the 
1857  pg/cm  titanium  load  is  that  the  K-shell  emission  becomes  very  model  sensitive  for  loads  that 
have  values  of  rj*  <  2.0.  A  value  of  rj*  near  unity  implies  that  there  is  barely  enough  energy  to  ionize 
into  and  radiate  from  the  K  shell,  and  its  value  is  1.5  for  this  titanium  load.  Titanium  experiments 
performed  on  the  Z  machine  also  exhibited  a  decline  in  K-shell  yield  in  this  rj*  regime.15 

Based  on  the  good  K-shell  yield  agreement  between  the  TCRE  model  and  the  “ on-the-spot” 
model  for  the  1857  pg/cm  titanium  load  and  the  disagreement  between  the  “ on-the-spot ’  and 
GCRE  model,  it  is  reasonable  to  speculate  that  the  major  deficiency  of  the  TCRE  modeling  of  these 
load  is  the  use  of  the  “ on-the-spot ”  approximation.  Therefore,  in  the  future,  it  may  be  worthwhile 
to  improve  the  TCRE  model  by  using  an  escape  probability  that  is  based  on  a  more  accurate  transit 
probability  than  that  calculated  using  the  “ on-the-spot ”  approximation. 

The  ion  density  and  electron  temperature  profiles  produced  by  the  TCRE,  “ on-the-spot' ,  and 
GCRE  models  were  generally  in  good  agreement  for  the  four  tested  loads.  The  greatest  deviation 
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occurred  for  the  500  pg/cm  argon  load.  The  profiles  for  this  load  are  shown  in  Figures  4  and  5 
for  the  ion  density  and  electron  temperature,  respectively.  The  time  profile  of  the  K-shell  radiative 
power  is  also  shown  for  this  load  in  Fig.  6.  Even  the  profiles  for  the  1857  pg/cm  titanium  load, 
which  had  the  largest  model  disagreement  in  terms  of  predicted  K-shell  yield,  showed  generally 
good  profile  agreement.  This  is  illustrated  in  Figures  7-10,  respectively,  for  the  ion  density,  electron 
temperature,  total  radiative  power,  and  K-shell  power  profiles.  For  clarity,  the  “on-the-spot’  results 
are  omitted  from  these  figures.  They  lie  between  the  TCRE  and  GCRE  results.  In  all  cases  there 
were  significant  spatial  and  temporal  differences  between  the  profiles  generated  under  optically 
thin  conditions  and  the  profiles  generated  using  any  of  the  other  models. 

Table  II  shows  a  typical  example  of  the  computer  times  that  the  GCRE,  TCRE,  and  “on-the- 
spot”  models  require  to  perform  a  calculation.  Specifically  the  timing  results  for  the  1857  pg/cm 
titanium  load  calculation  are  displayed.  The  time  spent  transporting  the  radiation  dominates  the 
total  calculation  time  for  the  GCRE  and  “on-the-spot’  models,  whereas,  the  radiation  computation 
time  is  negligible  for  the  TCRE  model.  The  radiation  calculation  was  at  least  1000  times  faster  for 
the  TCRE  model  than  for  either  the  GCRE  or  “ on-the-spot ’  models.  Note,  these  calculations  were 
performed  using  a  500  MHZ,  Pentium  III  personal  computer. 

In  the  2D  calculations  discussed  below,  the  radiation  computation  time  for  the  TCRE  model 
is  about  one  third  of  the  total  computation  time.  The  reason  this  time  is  relatively  longer  in  2D  than 
ID  is  because  of- the  2D  ray  trace.  For  a  fixed  number  of  rays  the  computation  time  associated  with 
the  2D  ray  trace  scales  approximately  with  the  square  of  the  number  of  zones.  Furthermore,  the 
number  of  rays  required  to  adequately  resolve  the  radiation  field  is  much  greater  in  2D  ( typically 
50  to  100  rays  are  needed). 

IV.  Two-Dimensional  TCRE  Calculations 

The  main  purpose  of  this  section  is  to  demonstrate  that  the  TCRE  model  can  be  readily  adapted 
to  an  existing  2D  code,  namely  PRISM,  and  to  present  some  2D  results.  However,  we  also  consider 
some  modeling  issues  in  ID  and  2D  that  can  have  as  profound  an  effect  on  results  as  the  choice 
of  radiation  model.  A  proven  method  of  increasing  understanding  and  improving  models  is  to 
compare  experimental  data  with  model  calculations.  By  improving  the  radiation  models  used  in 
the  calculations  one  can  better  focus  on  other  physics  and  numerical  issues  to  resolve  differences 
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between  experimental  and  theoretical  results. 

A  comparison  is  made  between  ID  and  2D  model  results  for  K-shell  radiative  yields  obtained 
from  argon  gas  puff  calculations  that  were  intended  to  simulate  experiments  performed  on  the  Z 
machine  at  Sandia  National  Laboratories.  These  gas  puff  calculations  are  for  2  cm  length,  7  cm 
diameter,  double  puffs  of  500,1000,  1500,  2000,  and  2500  pg/cm.  The  mass  distribution  of  gas 
used  in  these  calculations  is  taken  from  the  measured  gas  distribution  obtained  from  interferometry 
measurements  performed  at  the  Naval  Research  Laboratory,16  see  Figure  11.  A  circuit  model  of 
the  Z  machine  is  used  to  drive  the  theoretical  implosions. 

The  mass  distribution  is  initially  uniform  along  the  z  direction  in  the  2D  calculations;  however, 
when  the  outer  radius  of  the  pinch  has  imploded  to  a  radius  of  less  than  3  cm  a  density  and 
velocity  perturbation  is  applied.  The  initial  perturbation  wavelength  A  is  chosen  as  0.5  cm,  which 
is  a  typical  disruptive  wavelength.  It  is  a  non-compressible  perturbation  characterized  by:  (1) 
p'  =  pQ  x  (0.25)  x  cos{2n  z  /  X)\  where  the  mass  density  p  equals  p0  (initial  mass  density)  +p'  and 
(2)  vl  =  —u0  x  (0.25)  x  cos(2ttz/\)-,  where  the  radial  velocity  u  equals  u0  (initial  radial  velocity) 
+  u' .  When  the  perturbation  is  turned  on  it  is  limited  to  zones  near  the  outer  plasma  radius  because 
it  is  only  applied  when  |u|  >  0.1  x  \umax\.  Note,  the  choice  of  perturbation  level  and  tum-on  time 
is  somewhat  arbitrary,  which  is  in  the  spirit  of  this  preliminary  investigation  to  look  for  general 
trends  made  apparent  by  the  2D  modeling  of  Rayleigh-Taylor  instability  growth.  Figure  12  shows 
the  2D  radiative  power  density  profile  for  a  2000  pg / cm  load  at  30  ns  after  the  perturbation  tum-on 
time.  Notice,  that  this  figure  shows  that  the  radiative  power  is  larger  at  the  axial  boundaries,  which 
is  a  consequence  of  their  higher  photon  escape  probabilities. 

The  K-shell  yield  comparison  is  made  for  the  following  models:  (1)  ID  Lagrangian  MHD 
with  GCRE  radiation  transport,  (2)  ID  Lagrangian  MHD  with  TORE  radiation  transport,  (3)  ID 
Lagrangian  MHD  with  optically  thin  radiation,  (4)  ID  Eulerian  MHD  (PRISM  code)  with  TCRE 
transport  and  (5)  2D  Eulerian  MHD  (PRISM  code)  with  TCRE  transport.  The  ID  Lagrangian 
calculations  employ  80  radial  zones,  the  ID  PRISM  calculation  employs  200  radial  zones  and  the 
2D  PRISM  calculation  uses  200  radial  zones  and  40  axial  zones.  The  only  substantial  difference 
between  applying  the  TCRE  model  in  2D  as  opposed  to  ID  is  that  the  ray  trace  method  used 
to  calculate  the  “ on-the-spot' ’  escape  probabilities  is  more  complicated.  The  2D  method  used 
here  was  developed  by  Chong17.  Nearly  the  same  K-shell  yield  and  total  yield  results  (within  10 
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percent)  are  attained  with  the  ID  Eulerian  code  when  it  is  run  with  either  the  same  ray  trace  scheme 
utilized  in  the  ID  Lagrangian  calculations  or  with  the  2D  ray  trace  method.  The  K-shell  yields  as 
a  function  of  argon  mass  load  are  plotted  for  these  five  models  in  Figure  13. 

The  ID  Lagrangian  GCRE  K-shell  yield  results  curve  (1)  of  Fig.  13  is  nearly  identical  to  the 
ID  Lagrangian  TCRE  curve  (2).  This  agreement  and  the  fact  that  the  n*  and  Te  stagnation  profiles 
produced  by  these  two  models  are  also  in  good  agreement  (not  shown),  provide  strong  evidence 
that  the  TCRE  model  is  modeling  opacity  and  optical  pumping  effects  in  accord  with  the  GCRE 
model.  Therefore  the  TCRE  model  should  also  be  suitable  for  modeling  these  effects  in  2D.  The 
substantial  differences  between  the  GCRE  curve  (1)  results  and  the  optically  thin  results  of  curve 

(3)  again  illustrate  the  importance  of  including  detailed  and  accurate  radiation  models  in  Z-pinch 

calculations. 

Comparing  the  ID  Eulerian  results  of  curve  (4)  with  the  ID  Lagrangian  results  of  curve 
(2)  reveals  that  there  are  large  inherent  numerical  differences  when  implosion  calculations  are 
performed  in  a  Lagrangian  versus  an  Eulerian  reference  frame.  In  all  fairness  to  the  Eulerian 
model,  there  was  probably  inadequate  resolution  in  the  radial  direction  to  converge  its  results. 
Nevertheless,  the  point  needs  to  be  made  that  model  characteristics  other  than  the  choice  of 
radiation  treatment  can  affect  results.  In  particular,  model  reference  frame  and  zoning  can  have  a 
strong  influence.  Minimizing  some  of  the  model  dependencies  by  benchmarking  a  calculation  to  a 
specific  experiment  and  then  analyzing  the  rest  of  the  experiments  in  the  senes  under  the  conditions 
of  the  benchmark  is  a  useful  technique  for  studying  ID  and  2D  experimental  trends.13-18’19  This 
benchmarking  has  involved  alterations  in  transport  coefficients,  perturbation  levels,  perturbation 
tum-on  times,  and  radiation  dilution  factors. 

The  2D  Eulerian  results  of  curve  (5)  when  compared  with  the  ID  Eulerian  results  of  curve 

(4)  indicate  that  the  presence  and  growth  of  R-T  instabilities  in  the  experiments  has  the  potential 
to  somewhat  degrade  the  K-shell  yields.  However,  in  this  preliminary  investigation,  a  comparison 
of  these  curves  with  the  ID  Lagrangian  results  of  curve  (2)  again  emphasizes  that  the  numerical 
algorithm  and  zoning  can  play  at  least  as  large  a  role  as  R-T  instabilities  in  affecting  calculated 

K-shell  yields. 


V.  Summary 
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A  computationally  efficient,  multi-dimensional,  transport  method  suitable  for  modeling  Z- 
pinch  plasmas  that  are  in  collisional  radiative  equilibrium  has  been  developed  and  demonstrated 
in  ID  and  2D  calculations.  This  method  is  superior  to  diffusion  methods  because:  (1)  the  source 
function  is  based  on  a  self-consistent  non-LTE  calculation,  not  an  ad  hoc  assumption  that  it  is 
Planckian,  (2)  the  optically  thin  radiation  is  allowed  to  escape,  not  subjected  to  an  inappropriate 
diffusion  equation,  and  (3)  accurate  line-by-line  escape  probabilities  are  used,  not  diluted  as  in 
diffusion,  by  frequency  groups  that  mix  line  and  continuum  opacities. 

The  basis  of  this  TCRE  model  is  an  approximation  that  states:  a  localized  region  of  plasma 
with  ion  density  n*,  electron  internal  energy  E,  and  escape  probability  V  for  a  given  line  photon 
will  have  the  same  line  emission  coefficient  and  absorption  coefficient  as  that  given  by  a  uniform 
plasma  having  identical  n*,  E,  and  V.  This  approximation  allows  one  to  utilize  existing  ID  escape 
probability  based  results  for  the  solution  of  the  radiative  transfer  equation  to  approximate  a  solution 
in  two  or  three  dimensions.  The  TCRE  model  is  very  efficient  because  the  ID  transport  results  that 
it  is  based  upon  are  tabulated.  This  method  is  at  least  a  thousand  times  faster  than  the  time  required 
for  the  full  transport  calculation. 

The  current  implementation  of  the  TCRE  model  is  in  reasonably  good  agreement  with  the 
GCRE  calculation  over  a  wide  range  of  Z-pinch  titanium  and  argon  plasma  conditions,  even  when 
sharp  gradients  in  temperature  and  density  are  present.  We  have  not  encountered  conditions  where 
the  TCRE  approach  gives  poor  results.  However,  there  is  missing  physics  in  the  TCRE  model  which 
must  be  added  in  the  future  to  guarantee  its  continued  robustness.  For  instance,  continuum  opacity 
is  not  included  in  the  calculation  of  the  escape  probabilities,  nor  can  time-dependent  ionization 
dynamics  be  treated.  As  Z  -pinch  drivers  become  more  energetic  and  powerful,  simulation  of  their 
plasmas  will  increasingly  require  that  these  effects  be  included  in  any  realistic  modeling.  Another 
potential  weakness  of  TCRE  is  that  the  present  calculations  have  shown  that  for  plasmas  just  barely 
energetic  enough  to  ionize  into  the  K  shell,  the  “ on-the-spoE  approximation  begins  to  break  down. 
We  expect  to  address  these  issues  as  further  development  of  the  TCRE  approach  proceeds. 

It  must  be  noted  that  developing  a  reliable  link  between  ID  Lagrangian  and  2D  Eulerian  via  ID 
Eulerian  descriptions  of  a  radiating  plasma  is  not  simply  a  technical  issue  of  numerical  simulations. 
Rather,  it  has  a  direct  practical  significance  for  the  development  of  plasma  radiation  sources  (PRS). 
The  full  transport  (GCRE)  ID  Lagrangian  calculations  describe,  to  the  best  of  our  knowledge,  the 
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radiation  performance  of  PRS  under  the  "best  case"  conditions  where  instability  of  the  plasma  flow 
is  completely  suppressed.  Such  modeling  is  fully  adequate  for  studying  the  effect  of  radial  load 
structures  on  the  K-shell  radiation  yield,  e.g.,  to  evaluate  the  contribution  of  this  effect  apart  from 
the  axial  dependence  of  the  radiated  power  especially  in  gas-puff  loads.  In  two  dimensions,  the 
effects  of  structure  and  stability  become  inseparable.  A  2D  Eulerian  RMHD  simulation  allows 
inclusion  of  both  effects,  but  without  accurate  radiation  transport  it  lacks  a  reliable  reference  point 
with  which  to  evaluate  the  degree  of  hydrodynamic  modeling  inaccuracy.  One  reference  point  is 
provided  by  the  GCRE  Lagrangian  ID  model.  Combined  use  of  both  models  will  allow  us  to  better 
address  the  major  issues  of  PRS  load  design  and  optimization.  For  instance,  if  it  is  established  that 
under  certain  load  conditions  the  K-shell  yield  is  sensitive  mainly  to  the  radial  density/composition 
structure  of  the  load,  then  the  instability  is  tolerable,  and  we  should  optimize  the  radial/composition 
structure,  e.g.,  through  gas  puff  nozzle  design.  On  the  other  hand,  if  the  key  factor  turns  out  to  be 
the  RT  instability  of  implosion,  then  this  is  the  issue  to  be  addressed  by  structuring  the  load  for  the 
specific  purpose  of  mitigating  this  instability,  e.g.,  by  adding  a  weak  axial  magnetic  field. 
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Fig.  1.  Schematic  diagram  of  the  Tabulated  CRE  transport  model. 

Fig.  2.  Electron  temperature  (Te  (eV))  and  effective  charge  (Z)  as  a  function  of  electron  internal 
energy  and  plasma  thickness  A r.  The  argon  ion  density  n*  is  7.0  x  1017cm-3. 

Fig.  3.  Electron  temperature  (Te  (eV))  and  effective  charge  (Z)  as  a  function  of  electron  internal 
energy  and  plasma  thickness  A r.  The  argon  ion  density  n*  is  7.64  x  10 19cm-3. 

Fig.  4.  Ion  density  (nj  radial  profile  at  the  time  of  peak  argon  K-shell  emission  for  a  variety  of 
radiation  transport  models.  Results  are  for  a  500  f. ig/cm  argon  load. 

Fig.  5.  Electron  temperature  (Te)  radial  profile  at  the  time  of  peak  argon  K-shell  emission  for  a 
variety  of  radiation  transport  models.  Results  are  for  a  500  \igjcm  argon  load. 

Fig.  6.  Time  evolution  of  the  K-shell  power  for  a  variety  of  radiation  transport  models.  Results 
are  for  a  500  ng/cm  argon  load. 

Fig.  7.  Ion  density  («*)  radial  profile  at  the  time  of  peak  titanium  K-shell  emission  for  a  variety  of 
radiation  transport  models.  Results  are  for  a  1857  fj-g/cm  titanium  load. 

Fig.  8.  Electron  temperature  (Te)  radial  profile  at  the  time  of  peak  titanium  K-shell  emission  for  a 
variety  of  radiation  transport  models.  Results  are  for  a  1857  fig /cm  titanium  load. 

Fig.  9.  Time  evolution  of  the  total  radiated  power  for  a  variety  of  radiation  transport  models. 
Results  are  for  a  1857  g,g/cm  titanium  load. 

Fig.  10.  Time  evolution  of  the  K-shell  radiated  power  for  a  variety  of  radiation  transport  models. 
Results  are  for  a  1857  \xg/ cm  titanium  load. 

Fig.  11.  Initial  ion  density  distribution  used  in  2D  PRISM  double  puff  calculations  to  preliminarily 
investigate  argon  experiments  performed  on  the  Z  machine  at  Sandia  National  Laboratories. 

Fig.  12.  Spatial  profile  of  the  total  radiative  power  30  ns  after  the  perturbation  is  turned  on.  Result 
is  for  a  2000  fig /cm  argon  double  puff  load  driven  by  a  Z  machine  circuit  model. 


48 


Fig.  13.  Comparison  of  K-shell  yields  as  calculated  from  ID  and  2D  models  of  argon  gas  puff 
experiments  performed  on  the  Z  machine.  Calculations  are  performed  for  mass  loads  of  500, 1000, 
1500,  2000,  and  2500  \igjcm. 
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Table  I.  Radiation  transport  model  comparison  of  total  and  K-shell 
radiative  yields  in  (kJ/cm)  for  different  Z-pinch  loads.  The  K-shell 
radiative  yields  are  in  parenthesis. 


Radiation 

Transport 

Model 

Argon 

150 

jLig/cm 

Argon 
500 
jLig/c  m 

Titanium 

1125 

JLig/cm 

Titanium 

1857 

jug/cm 

TCRE 

12  (3) 

87  (28) 

262  (70) 

366  (42) 

On-the-Spot 

10  (3) 

81  (28) 

258  (73) 

370  (39) 

GCRE 

10  (3) 

74  (23) 

263  (69) 

395  (26) 

Optically  Thin 

23  (3) 

100  (30) 

400  (20) 

600  (4) 
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Table  II.  Radiation  transport  model  comparison  of  computer  timings 
for  the  ID  1857  mg/cm  titanium  wire  array  calculation. 


Radiation 

Transport 

Model 

Total  calculation 
time  (cpu  -  seconds) 

Radiation  transport 
calculation  time 
(cpu  -  seconds) 

GCRE 

29270 

26286 

"On-the-spot" 

19280 

16296 

TORE 

2984 

11 
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1  -  Optically  thin 

2  -  TORE 

3  -  "On-the-spot" 

4  -  GCRE 
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IV.  Wire  Load  Dynamics  and  Modeling 


for  DTRA  Plasma  Radiation  Sources 


Wire  loads  have  shown  a  productive  and  reliable  path  to  record  high  PRS 
radiative  yields  and  powers.  The  technology  offers  a  heretofore  unknown 
level  of  control  over  pinch  impedance  history,  load  material  variations,  and 
assembly  processes.  As  an  element  of  the  DTRA  simulator  program,  the 
research  and  development  effort  focused  on  a  fundamental,  predictive  and 
versatile  understanding  of  these  loads  is  described  here. 

An  early  insight  into  the  operation  of  a  wire  load,  in  either  a  single  or 
nested  array,  was  provided  by  the  advent  of  the  Wire  Dynamic  Model  (WDM) 
wherein  the  wire  elements  are  treated  as  distinct  circuit  elements.  The  model 
has  successfully  described  mass  and  current  interchange  in  nested  array  loads 
over  a  current  domain  of  0.1  to  20  MA  and  the  mass  ranges  that  go  with 
that  for  implosion  times  from  100  to  300ns.  The  behavior  of  bias  induc¬ 
tance  and  wire  resistance  has  been  investigated  experimentally  and  model 
generalizations  to  treat  those  variations  are  discussed  here. 

On  the  other  hand  a  closer  examination  of  wire  implosion  dynamics  and 
initiation  processes  reveals  a  rich  set  of  phenomena  that  are  not  so  well 
understood.  One  such  feature  is  the  axial  precursor  of  a  wire  array  which 
forms  on  axis  as  array  material  is  sloughed  off  and  forced  inward  by  magnetic 
pressure.  The  precursor  is  remarkably  stable  and  settles  into  a  steady  state 
radius  which  is  apparently  characteristic  of  the  the  material  and  its  radiation 
properties.  A  second  such  feature  is  the  generation  of  a  wire  corona  plasma 
which  feeds  the  precursor  and  controls  the  inertia  of  the  parent  wire  array 
as  the  implosion  progresses.  A  third  related  feature  is  the  process  by  which 
a  wire  array  “closes”  electrically  and  the  relative  fraction  of  current  that  is 
carried  in  the  denser  portions  of  the  wire  plasma  at  the  time  this  closure 
obtains. 

A  fair  amount  of  research  to  date  appears  to  indicate  that  tenuous  plasma 
dynamics  is  a  noticeable  to  dominant  contributor  to  the  wire  array  implosion 
physics  over  the  operational  range  of  wire  masses,  interwire  gaps,  and  deliv- 
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ered  current  per  wire.  It  is  also  clear  that  the  quality  of  equation  of  state 
and  resistivity  data  plays  a  major  role  in  the  quality  of  theoretical  prediction 
of  the  wire  array  physics. 

Treating  the  tenuous  plasma  stresses  conventional  MHD  codes  nearly  to 
(or  argueably  past)  the  point  of  validity  as  such  thin  material  either  streams 
around  the  interior  wires  of  a  nested  array  or  builds  a  precursor  column. 
The  corona  plasma  experiences  a  wide  dynamic  range  in  ionization,  perhaps 
radiatively  driven,  and  conductivity  as  well.  A  final  topic  of  research  dis¬ 
cussed  here,  the  evolution  of  the  WDM  toward  a  true  fluid  particle  model, 
is  motivated  primarily  by  the  requirement  to  understand  and  to  control  (or 
perhaps  even  exploit)  these  load  details  for  the  optimization  of  the  radiation 
product. 

While  the  results  of  this  effort  are  somewhat  formal  in  nature,  the  re¬ 
quired  elements  for  a  complete  model  are  now  in  hand  and  some  important 
properties  of  the  theory  are  demonstrated.  At  present  the  model  can  be  char¬ 
acterized  as  a  Drift  Kinetic  Fluid  Particle  Hydrodynamics  (DKFPH)  which 
will  span  densities  characteristic  of  dense  pinch  regions  and  corona  with  the 
same  algorithm.  The  new  fluid  elements  can  be  equipped  with  any  desired 
equation  of  state,  transport  coefficients,  and  modern  radiation  treatments. 
A  systematic  particle  management  scheme  is  derived  from  straightforward 
thermodynamic  considerations. 

1  Applications  Improving  the  WDM 

The  Wire  Dynamic  Model  is  a  special  case  of  a  full  fluid  MHD  picture  in 
that  (i)  the  wires  are  represented  only  as  particles,  while  (ii)  the  dissipative 
or  collisional  effects  in  these  wires  are  treated  in  a  spatially  averaged  way,  and 
(iii)  the  external  forces  are  developed  by  an  inductive  and  resistive  lumped 
circuit  description  that  presumes  a  known  spatial  concentration  of  electrical 
conductivity.  These  specializations  are  mostly  of  historical  significance  and 
are  sequentially  being  replaced  by  more  complete  physics.  However,  the  sys¬ 
tematic  perturbation  about  such  valuable  results  has  enabled  the  retention  of 
an  application  ready  version  while  the  better  implementations  are  developed. 

In  the  following  sections  are  discussed  several  new  results  which  were 
added  to  the  WDM  during  the  program  year,  some  are  primarily  formal 
while  others  are  focused  on  particular  experiments. 
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1.1  Energy  Conservation  and  the  Circuit  Model 

A  definitive  treatment  of  Poynting’s  theorem  for  the  wire  ensemble  was  not 
available  once  sticky,  interacting  wires  were  brought  into  play.  The  changes 
wrought  by  a  better  treatment  are  reflected  in  a  more  robust  conservation  of 
energy  and  the  ability  to  capture  the  trade  between  inductive  and  resistive 
current  division  that  is  being  observed  in  the  laboratory. 

For  the  WDM  viewed  as  a  “Biot-Savart  gas”  the  direct  interaction  among 
an  ensemble  of  current  elements  is  equivalent  to  a  potential  energy  defined 
by  the  array  inductance  and  the  velocity  dependent  drag  forces  axe  derivable 
from  a  Rayleigh  dissipation  function.  Within  the  ensemble  are  symmetry 
groups  of  multiplicity  S,  e.g. 


where  on  is  any  S- wise  invariant  current  fraction  contained  in  any  particu¬ 
lar  wire  path  represented  by  a  series  resistance  and  inductance.  Once  the 
trajectories  of  an  ensemble  of  wire  elements  axe  followed, 


IMX-is  =  ( 


27T 


-)C 


sp  K*ij  sp /vu 


-( «r>) 

determines  the  motion  and  the  time  advanced  impedances  L ij.  The  energy 
absorbed  and  the  work  expended  against  the  drag  forces  can  be  viewed  as 
deposited  in  two  reservoirs,  the  inductive  energy  £$  =  l/2I2(t)aiLi:,o;j  and 
the  dissipation  function, 

Tij{t  -  t0)  =  1/2  /‘r(^)(xy)2dt  , 

J  to 

with  Uij  a  radial  separation  and  aw  a  wire  size. 

By  construction  then  the  sum  of  all  on  is  one.  If  we  examine  the  voltage 
Ve(t)  impressed  at  the  entrance  or  feed  of  a  wire  array  cage,  then  from  the 
following  figure,  the  resolution  of  the  set  on  and  the  total  current  1  admitted 
to  the  array  are  seen  to  be  two  distinct  problems.  From  the  view  of  the 
external  generator  Vg  (t) ,  the  array  is  a  two  terminal  device  and  the  energy 
into  it  must  be  the  same  for  all  sets  a*  that  show  the  same  net  impedance. 
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WDM  circuit  elements  can  decompose  into  symmetry  groups, 
with  each  element  representing  a  resistance  and  inductance. 


Maintaining  equal  voltage  at  each  circuit  element,  the  equations  fixing 
the  current  fractions  a*  for  any  given  Ve(t)  =  Vg(t)  —  ZgenX  and  X  can  be 
written  as, 

and  then  integrated  over  the  interval  [<5  =  t>  —  t<]  to  eliminate  the  time 
derivative.  Denote  by  <5$e  the  change  in  flux  over  this  interval,  the  matrix 
relation 


6$e 

j  ”1”  hi,*  *  Ctj 

emerges,  with  Iy  the  diagonal  identity  matrix,  U*  a  unit  vector  over  the  local 
group  N.  The  common  LHS  to  all  these  equations  can  be  expressed  as  an 
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effective  impedance  at  the  entry  point,  viz.  5Ze  =  ^  for  any  impressed 
voltage  and  current.  One  may  also  write  Le  =  SZe  to  put  the  linear  system 
in  units  of  inductance.  Normalized  by  the  port  inductance  Le  the  equations 
for  the  current  fractions  are  simply, 

LeUi  =  (SRjlij  -t-  L ij)  •  ocj  , 

and  since  the  fractions  must  sum  to  one,  the  port  inductance  or  impedance 
may  be  calculated  independently  of  the  particular  value  of  voltage  or  current 
in  the  array  of  conduction  paths.  It  is  easy  to  show  that  in  fact 

Lg1  =  Uj  •  (6Rjlij  +  L ij)  ■  Uj  . 

In  the  limit  of  vanishing  inductance  it  is  easy  to  see  that  the  oci  =  Z\\J  Ri{t ) 
and  that  consequently  they  add  up  to  unity. 

In  the  general  case  the  diagonal  nature  of  the  series  resistances  breaks 
the  problem  of  energy  transfer  into  a  balance  between  the  power  building  the 
electromagnetic  energy  density  in  the  gap  and  the  power  brought  forward  of 
the  generator  impedance  and  series  wire  resistances.  On  the  left  hand  side 
is  the  Poynting  flow  feeding  the  wire  load,  with  power  removed  because  the 
series  resistors  kill  some  magnetic  flux, 

I  X  (ve(t)  —  TZiRi{t)af^  —  —j^(£field  d"  £mech)  > 

which  is  balanced  on  the  right  by  the  total  time  derivative  of  field  and  me¬ 
chanical  energy.  With  the  dissipation  function  defined  above,  we  can  inte¬ 
grate  over  the  previously  used  time  interval  5  to  obtain  the  change  in  energy 

AS  =  (: Ts  -  T)  -  S ^(S)  +  1/2 (LaS%  -  LJ?)  , 

with  5  indicating  the  foward  time  level  for  kinetic  energy,  T  =  1/2 
and  magnetic  energy,  1/2  LJ?  =  1/2  ^[EyOjLy^].  Of  course  here  the  on 
are  normalized  to  contain  all  the  current  in  the  ith  wire  of  all  the  symmetry 
groups. 

Time  integrating  the  delivered  power  requires  two  special  functions,  Ci(o) 
and  C2(5),  easily  developed  from  the  “L/R”  solution  waveforms  characteriz¬ 
ing  the  feed.  The  result  is  dependent  on  the  time  interval  only  through  the 
two  coefficients  and  becomes 

Vs  =  Cx(S)Vgl  -  C2(6)l\Zgen  +  Si  Riot)  ■ 
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Combining  this  result  with  the  expression  for  the  change  in  field  and  mechan¬ 
ical  energy  leads  to  an  algebraic  constraint  for  the  time  advanced  current  Z& 
as  the  solution  of  a  quadratic  equation,  viz. 

V5  =  C Ts  -  T)  -  E ijFijiS)  +  l/2(LaSl$  -  Lal 2)  . 

The  only  required  input  is  a  voltage  waveform  Vg(t).  When  carefully 
implemented  this  algorithm  allows  conservation  of  energy  to  essentially  ma¬ 
chine  precision  with  very  tight  limits  placed  on  the  changes  in  wire  velocity 
and  position.  The  trade  between  energy  convergence  and  runtime  as  such 
limits  are  relaxed  is  also  very  favorable.  A  relaxation  scheme  to  converge  the 
particle  positions,  velocities,  and  currents  together  is  completely  successful. 

The  accurate  exchange  of  field  and  mechanical  energy  is  a  cornerstone  of 
any  further  enhancements  of  the  WDM  that  evolve  it  to  a  full  fluid  particle 
description.  It  has  also  made  possible  the  accurate  assessment  of  how  impor¬ 
tant  flux  compression  events  are  to  the  nested  wire  loads  from  the  viewpoint 
of  a  filamentary  current  model. 

1.2  Imperial  College  Experiments  and  the  WDM 

By  way  of  review,  the  initial  successes  of  nested  arrays  motivated  an  attempt 
at  Imperial  College  to  understand  the  transfer  mechanisms  for  current  and 
momentum  from  an  experimental  perspective. 

The  nested  array  geometry  was  suggested  as  a  means  of  obtaining  a  fast 
rising  current  on  the  interior  (or  target)  array  in  J.  Davis,  et.  al.,  Appl. 
Phys.  Lett.  70(2),  13  Jan  97.  The  nested  geometry  was  found  to  enhance 
x-ray  power  on  the  “Z”  facility  at  20  MA.  The  best  yields  and  powers  were 
seen  with  a  target  array  lighter  than  the  driver  in  C.  Deeney,  et.  al.,  PRL 
81(4883),  30  Nov  98.  The  Wire  Dynamic  Model  showed  that  observed  pinch 
dynamics  on  “Z”  required  mass  interchange  and  current  switching  for  correct 
modeling  c.f.  R.  E.  Terry,  et.  al.,  PRL  83(4305),  22  Nov  99. 

Operating  at  a  lower  current  (1  M A),  detailed  observations  of  nested  array 
dynamics  on  MAGPIE  show  the  mass  interpenetration  and  show  a  shift  in 
assembly  time  due  to  added  inductance  between  the  drive  and  target  arrays, 
c.f.  S.  V.  Lebedev,  et.  al.,  PRL  84(1708),  21  Feb  00.  All  lower  current, 
lower  number  array  implosions  have  been  characterized  by  the  formation  of 
early  plasma  on  axis.  The  implosion  dynamics  are  strongly  modified  by  this 
behavior  —  S.  V.  Lebedev,  et.  al.,  PRL  81(4152),  9  Nov  98,  J.  P.  Chittenden, 
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et.  al.,  PRL  83(100),  5  July  99.  Here  we  seek  to  examine  these  observations 
against  the  numerical  model. 


1.2.1  Comparing  the  Model  Dynamics  to  Experiment 

The  experimental  geometry  is  shown  in  this  figure,  with  the  offset  inductance 
indicated. 


Magpie  Front  End  and  Equivalent  Circuit 

Electrical  and  mechanical  data  due  to  S.  V.  Lebedev 


The  relevant  parameters  for  the  simulation  of  implosion  events  in  this 
device  were  taken  as:  i.  array  diameter  of  8  mm  Target  and  16  mm  Drive, 
ii.  wire  diameter  of  1.50-10-3  cm,  iii.  wire  length  of  23  mm  +  offset  of  0  ^ 
33  mm  in  Target  array,  and  iv.  an  array  mass  (D  or  T)  of  0.1748  mg. 

The  WDM  model  was  run  against  a  series  of  test  problems  with  varying 
offset  length  because  the  amount  of  conductor  remaining  in  the  circuit  when 
the  target  wires  break  is  unknown.  The  high  and  low  inductance  cases  shown 
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below  axe  on  the  extreme  range  of  the  study. 


Early  Mass  Shedding  to  a  Precursor  is  Clearly  Important  The 
WDM  at  fixed  element  mass  fails  to  get  the  right  strike  time  or  trajectory 
shape  before  the  current  transfer.  On  the  other  hand,  the  mass  shedding  isn’t 
important  after  the  strike  —  we  get  a  good  trajectory  calculation  without  it. 
The  current  transfer  to  the  target  array  must  therefore  be  much  faster  than 
the  risetime  in  the  main  circuit. 

High  Inductance  Nested  Array  Implosion 


Most  of  the  Offset  Inductance  Remains  in  the  Circuit  When 
the  offset  inductance  is  added,  the  best  results  are  obtained  when  a  full  33 
mm  of  length  is  used.  Moreover,  this  inductance  appears  to  remain  in  the 
circuit  after  the  strike  —  changes  in  this  value  for  the  offset  length  take  the 
calculation  away  from  the  data. 
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Current  [MA] 


Flux  Compression  is  not  important  in  MAGPIE  True  flux  com¬ 
pression  is  signaled  by  a  current  partition  fraction  aj  <  0.  In  modeling 
regular  experimental  cases  for  MAGPIE,  such  events  are  prolonged  and  no¬ 
ticeable  only  when  an  interior  wire  element  rebounds  from  the  origin  the 
“roach  motel”  effect.  We  can  make  the  effect  more  pronounced  during  the 
rundown  by  examining  azimuthally  modulated  cases  which  show  strong  fila- 
mentation  modes,  but  it  is  still  a  fairly  isolated  behavior. 


Low  Inductance  Nested  Array  Implosion 
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Direct  Momentum  Transfer  at  Strike  is  Low  In  the  experience 
to  date,  the  use  of  “sticky”  wires  will  cause  strong  divergence  of  our  calcu¬ 
lated  trajectories  from  the  observations.  Apparently  most  of  the  drive  array 
material  moves  through  to  the  origin  without  pushing  on  the  target  array 
directly. 

These  observations  are  distilled  into  the  best  fit  of  model  to  experiment, 
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Current  [MA] 


1.2.2  Conclusions 

Prom  the  WDM  analysis  of  this  experiment  several  observations  emerge. 

Mass  Interchange  and  Current  Switching  An  understanding  of 
these  effects  in  nested  array  dynamics  is  critical  to  a  correct  account  for 
array  currents  from  1  MA  to  over  20  MA. 

Flux  Compression  from  a  Rebounding  Drive  Array  The  “roach 
motel”  effect  on  rebounding  drive  plasma  once  it  gets  inside  the  target  array 
may  be  an  important  design  option  to  be  examined  by  interior  array  wire 
number  scans  on  bigger  machines. 


2  Precursor  Plasmas  in  Steady  State 

It  has  been  observed  [Lebedev,  et  al]  that  heavier,  sparse  wire  arrays  do 
not  simply  implode  under  the  influence  of  the  injected  current.  The  general 
behavior  is  the  creation  of  corona  plasma  around  a  core  of  wire  material,  with 
the  coronal  plasma  continually  streaming  into  the  center  of  the  array,  and 
impinging  on  the  interior  or  target  array.  This  core  and  corona  morphology 
leads  to  two  primary  consequences. 

First,  the  wire  cores  do  not  implode  following  the  dynamics  of  isolated 
wire  filaments.  The  mass  shed  to  the  corona  plasma  takes  on  the  momentum 
transferred  from  the  pulse  power  feed  and  leaves  the  core  to  be  smoothly 
boiled  away  for  a  majority  of  the  current  pulse.  Only  after  a  significant  frac¬ 
tion  of  the  original  mass  is  shed  does  the  set  of  wire  cores  begin  to  implode. 

Second,  the  coronal  plasma  forced  to  the  interior  generally  forms  a  very 
regular,  stagnant  cylinder  of  progressively  denser  plasma  that  radiates  away 
all  the  energy  brought  to  it  by  the  imploding  jet.  It  is  observed  that  the 
radius  of  the  column  depends  strongly  on  the  wire  material  used,  but  only 
weakly  on  time.  No  outwardly  moving  shock  front  is  observed,  but  rather  a 
progressively  denser  column. 

2.1  The  Central  Column 

The  very  regular  behavior  of  the  system  naturally  leads  to  the  problem  of  a 
radiative  calibration  of  the  precursor  plasma,  viz.  can  the  absolute  radiative 
output  be  uniquely  related  to  the  flux  contained  in  the  input  jet?  Can  the 
observed  variation  of  stagnation  radius  with  atomic  number  be  explained? 

As  a  starting  point,  we  will  consider  first  the  construction  of  a  simple 
steady  state  flow  model  which  describes  the  transport  of  mass  to  the  stagnant 
central  column. 

2.2  Geometry 

As  shown  in  the  figure,  the  flow  we  seek  to  describe  exists  only  on  an  annulus 
r,  <  r  <  r0.  At  the  outer  boundary  there  is  a  known  flow  speed  u0,  density 
no, "and  entrained  current  J0  which  is  generally  a  small  fraction  of  the  current 
in  the  drive  array  that  generates  this  inward  streaming  corona.  The  boundary 
mass  flux  T0  =  r0u0n0  =  T(C)  is  presumed  fixed  in  time  and  space  over  the 
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annular  domain,  thus  serving  to  define  the  interior  flow  profile  U(()  once  the 
interior  density  profile  N(()  is  known. 

Fig.  1  Steady  State  Adiabatic 


Flows  on  an  Annulus.  no  |o  po 


U(0  =1/CN(0 
.  •  r0=Ron0u0 

2.3  Model  Statement 

Mass  continuity  is  satisfied  by  construction,  viz. 

dtrMO  =  -C_1acr(C)  =  0  , 

with  £  =  r/ro.  Momentum  balance  is  to  be  established  from  density,  current, 
and  pressure  profiles.  The  density  form  factor  is  just  N(()  =  n(C)/n0  such 
that  N(l)  =  1.  The  current  is  derived  from  a  constant  current  density 
jz  =  jQ  .  T  =  (cr0Ezo)  •  T \  with  T  the  ionization  fraction  and  <To  the  dc 
conductivity  based  on  electron  neutral  collisions  at  small  T .  For  low  interior 
magnetic  field  and  flow  speeds,  the  transformation  to  the  moving  fluid  frame 
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is  a  perturbation,  as  is  the  spatial  dependence  of  Ez 0.  The  pressure  is  derived 
from  an  adiabatic  rule,  viz.  P(C)  =  PoN7(0-  These  assumptions  lead  to  a 
statement  of  momentum  balance  written  as 


-udru  =  p~ldcP  +  , 

47T  pQ 

with  Be  =  2/(C)/croC,  /(C)  =  W(C,C*).  and  /o  =  ^fco-Eo-  Since  the 
ambient  electric  field  is  assumed  constant  over  the  radial  domain  of  interest, 
the  current  density  thus  will  vary  in  accordance  with  the  field  in  the  radially 
moving  frame  that  is  seen  outside  of  £*,  viz. 


defines  the  current  density  form  factor.  Alternately,  one  may  cast  the  equa¬ 
tion  into  a  form  using  the  imbedded  current  form  factor  X,  that  we  defined 
above, 


Z(C,M=(C2-CV2^>J|y, 

Using  the  current  profile,  let  (B$  =  BqPI(C,  C»)  with  Bo  =  2/o/cro, 
and  set  <?s  =  Po/po-  Next,  writing  v2A0  =  B%/4np0,  and  substituting  u  = 
uo/CNiQ  one  finds  that  the  pressure  balance  relation  can  be  expressed  in 
terms  of  only  four  distinct  speeds,  or  equivalently  three  dimensionless  ratios 
—  a  =  cju o,  5  =  and  /?  =  vA0/u0.  In  particular, 


idrX 


~'=m9ifP+v^em' 

and  in  order  to  cast  the  relation  into  a  form  defining  C-N’(C)  rather  than  N(() 
only,  the  following  identity 


<9f  (C-7  exp(70(C)))  =  7C  7  exp(7^(C))  (d<G  -  ^  , 

is  useful.  Expanding  the  velocity  derivative  and  transforming  to  a  new  func¬ 
tion  £7(C)  =  ln[CN(C)]>  yields  a  first  order  system  of  ordinary  differential 
equations 


_a2TCe(7+1)g  +  egC7  01 

^7+1  _  (a027C^7+1i® 
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<9CJ  =  2  K  -  T6—^—r\  , 

C  V  exp  (Q)J 

where  we  have  collected  all  derivative  terms  to  the  left  hand  side.  These 
equations  are  now  to  be  solved  subject  to  the  boundary  conditions  Q((  — > 
1)  =  0  equivalent  to  N((  — ►  1)  =  1  since  exp[£(C)]  =  C/V(C),  and  J(C  — > 

Q  =  o. 

2.4  Model  Solution 

In  the  case  that  6  — »  0,  the  current  profile  can  be  explicitly  resolved  and 
only  a  single  ODE  remains.  Small  <5  is  an  appropriate  limit  when  /3  is  also 
small,  because  the  current  plays  a  minor  role  in  establishing  the  flow  profile. 
In  the  three  cases  show  below,  the  current  is  progressively  increased  and 
the  pressure  profile  'P(C)  tends  to  peak  toward  the  interior  as  a  measure  of 
the  stress  imposed  by  the  magnetic  field.  The  other  profiles  shown  (density 
AZ’(C),  flow  U(C),  and  G(Q)  are  as  defined  above. 
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When  the  ordering  a  >  f3  >  5  is  violated,  the  profiles  are  best  calculated 
by  the  full  pair  of  ODE’s  shown  above.  In  a  typical  case,  the  usual  condition 
of  constant  electric  field  in  the  lab  is  relaxed  and  a  parabolic  profile  rises  up 
from  a  selected  scale  radius  of  0.6. 


Here  the  current  density  J  is  extinguished  at  a  scale  radius  of  0.6,  while 
the  imbedded  electric  field  £  rises  to  the  value  required  to  ensure  that  no 
reversal  of  current  density  is  obtained  in  the  solution.  Once  such  a  balance 
is  struck  the  material  coasts  into  the  origin  with  a  very  flat  density  profile 
M  and  modest  upturn  of  velocity  U  as  it  nears  the  stagnation  zone. 

Using  a  decaying  electric  field  profile  serves  to  contain  the  current  to  the 
exterior  of  the  flow  in  keeping  with  observations  that  little  or  no  kink  mode 
instability  is  seen  in  the  stagnation  column. 
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2.5  Match  Conditions  at  C* 

Let  Qt  be  the  value  at  the  interior  boundary.  The  required  match  to  a 
stagnation  flow  region  there  is  determined  by  the  value  and  gradients  of  the 
fundamental  solution  curves. 

a  r  -o:27 Ce('y+1)g* 

°<,y*  -  ^7+i  _  (a^)27e(7+i )s.  ■ 


—  2£»  , 

and  the  implied  constructs  from  them:  M,  —  CT1^*)  U*  =  e  V*  = 
2c?elS'/Q,  U*  =  e~20'  +  2a2e(7-1)5* /C7-1  >  constituting  the  dimensionless 
form  factors  for  density,  velocity,  pressure,  and  specific  enthalpy  as  evaluated 
at  the  stagnation  boundary. 

At  the  stagnation  boundary  several  conditions  must  coexist.  First,  the 
flow  must  change  over  from  a  steady  flow  in  time  and  space  to  a  homoge¬ 
neous  compression  flow  which  collects  material  at  the  rate  determined  by 
the  steady  state  number  balance  condition.  Second,  the  mass  accretion  from 
the  thin  precursor  to  the  much  denser  central  column  must  conserve  radially 
directed  momentum  density  as  that  inelastic  collision  proceeds.  Finally  a 
power  balance  must  be  struck  between  the  radiative  losses  and  the  stagna¬ 
tion  enthalpy  flux  from  the  incoming  steady  state  flow.  Neither  the  first  nor 
second  condition  constrains  the  location  of  the  stagnation  radius,  but  it  is 
easy  to  see  that  the  third  determines  it,  and  hence  the  implied  column  size, 
uniquely. 

Turning  now  to  the  first  two  conditions,  when  the  interior  column  density 
profile  is  essentially  uniform  the  velocity  profile  will  be  linear  in  radius  and 
characterized  by  a  rate  parameter  that  itself  may  be  a  function  of  time,  viz. 
Kc  _  —aQ/2.  Now  the  density  will  build  according  to  the  simple  rule,  dtMc  = 
( uod/rQ)Afc .  As  fresh  material  is  added  in  an  inelastic  collision  process,  the 
radial  momentum  density  will  be  conserved  (provided  the  collision  is  confined 
to  a  small  radial  domain  and  no  external  forces  act).  Equating  the  interior 
momentum  density  to  the  column  momentum  density,  one  has  (A/'*+A/’c)£4  = 
AftU*,  which  leads  to  a  snowplow  constraint  on  the  recoil  speed  parameter, 


;  .  2  K 
a-UX  K+K  ’ 


with  the  density  evolving  according  to 
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dtNc=Ut~ 

Oo 


MMc 
A Cc  +  A/'c 


With  the  substitution  W  =  ln(Ac/A/"*),  the  density  is  to  be  calculated  as  a 
solution  to 


=  Z4 


2no  1 
C*ro  1  + 


with  the  boundary  condition  W(0)  =  0.  A  typical  solution  is  shown  in  the 
figure,  with  exponential  growth  giving  rise  to  a  quasi-linear  phase  after  a  few 
transit  times  r0/u0. 


Within  the  column  the  radiation  losses  will  scale  as  a  quadratic  to  lin¬ 
ear  function  of  (*Afc  depending  on  the  opacity,  and  the  spectrally  averaged 
emissivity  will  depend  most  sensitively  on  the  column  temperature  and  the 
atomic  species.  Hence  if  we  write  the  radiated  power  per  unit  length  in  terms 
of  the  form  factors  above,  VT  =  e(Tc,  A)(C*ew)2.  The  function  e(TC)  A)  is  also 
going  to  show  a  dependence  upon  the  detailed  column  density  profile  and  on 
the  rate  at  which  jet  energy  is  thermalized  (in  the  early  part  of  the  accretion), 
but  such  effects  are  less  important  than  the  scaling  with  atomic  species. 
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Since  the  column  will  fill  uniformly  with  time  as  the  precursor  pours  into 
the  center  the  radiative  power  might  be  expected  to  saturate  in  time  after  the 
material  becomes  more  optically  thick.  In  addition  at  later  times  that  optical 
thickness  will  make  the  column  more  of  a  surface  radiator  and  the  density 
dependence  will  cease.  However  these  limits  are  not  expected  in  the  case  of 
the  observations  due  to  Lebedev,  because  the  filling  rate  is  small  (r0/u0  may 
be  on  the  order  or  /as)  and  the  observed  radiation  is  not  that  of  a  blackbody. 

In  any  case  two  observations  are  now  crucial  to  a  description  of  the  stag¬ 
nant  precursor  column.  First,  every  relation  thus  far  derived  places  no  con¬ 
straint  whatever  on  the  parameter  («•  The  exterior  steady  state  flows  can 
show  a  wide  range  of  interior  end  points  as  a,  (3,  and  5  are  varied.  The 
collection  of  material  to  satisfy  number  balance  can  find  a  solution  for  any 
value  of  stagnation  radius.  Second,  nothing  in  the  steady  state  solution  is 
strongly  dependent  upon  the  atomic  species  involved.  Noticeable  sensitivity 
to  atomic  species  enters  only  through  the  radiation  loss,  e(Tc,  A)- 


Ke(Tc,  i4)(C*ew)2  =  e"25*  +  2e? e^l)9' / Q~l  , 

which  most  assuredly  does  constrain  £*.  The  scaling  factor  1Z  is  just  the 
ratio  of  the  radiation  power  at  AC  and  Tc  to  the  input  enthalpy  flux  rou0/io  — 
r0monoUo/2  before  any  radial  acceleration  of  the  fluid  occurs.  Generally  7Z  is 
a  large  number  and  so  the  solution  to  this  relationship  for  generally  tends 
to  occur  at  a  small  radius.  Moreover,  as  indicated  in  the  figure  above  that 
represents  a  typical  case,  larger  atomic  number  forces  a  solution  at  lower 
radius. 

In  the  cases  shown  above  the  upper  pair  is  from  a  set  of  profiles  with  a 
constant  electric  field  profile  and  no  6 .  The  lower  is  from  the  more  realistic 
choice  of  parabolic  electric  field  and  modest  5. 

3  Precursor  Investigations  with  Mach2 

The  active  version  used  here  was  Mach  2.5b  with  some  detailed  enhancements 
due  to  M.  Friese  and  M.  Douglas.  The  Lee-More-Desjarlais  resistivity  model 
for  A1  was  incorporated  with  a  set  of  SESAME  tables  build  from  data  due 
to  M.  Desjarlais. 

3.1  Problem  Statement 

Here  a  single  A1  wire  (at  2.0  cm  radius)  is  presented  with  a  scaled  Saturn  “log 
pulse”"  current  waveform  and  must  interact  with  a  return  current  conductor 
at  a  fixed  outer  radius  (2.5  cm).  The  lower  radius  of  the  problem  (1  cm) 
is  also  set  with  a  floating  conductor  boundary.  Symmetric  lateral  boundary 
conditions  are  imposed  appropriate  to  a  wire  count  of  24  or  so.  The  initial 
wire  is  set  to  100  fj,m  radius  to  allow  for  plenty  of  material  blowoff  into  a 
corona. 

3.2  Results 

In  keeping  with  observations  from  other  workers,  the  wire  tends  to  hold 
current  very  close  and  not  achieve  a  sufficiently  hot  corona  to  set  up  an 
inward  stream  of  material  to  feed  a  precursor.  The  magnetic  pressure  map 
shows  that  a  corona  would  indeed  drift  to  the  center  if  it  could  form. 
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LMD  RESISTIVITY 

Electric  field  structure 


I 

I 
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From  a  cold  start  the  Al 
wire  expands  with  the 
current  density  in  the  hot, 
dense  region  at  about  3.4 
eV.  The  exterior  is  about 
0.7  eV. 


4  Advances  in  Fluid  Particle  Methods 


The  general  success  of  the  wire  dynamic  model  (WDM)  leads  naturally  to¬ 
ward  a  whole  class  of  methods  based  upon  discrete  conduction  paths  and  fluid 
elements  that  can  more  easily  accomodate  the  phenomena  of  mass  interpen¬ 
etration,  current  exchange,  and  multiple  materials  for  radiation  sources.  It 
is  clear  that  azimuthal  structures  play  an  important  role  in  existing  load 
designs  and  may  in  fact  be  important  factors  in  new  load  designs  that  seek 
even  better  yields  and  more  controlled  performance  from  the  plasma  radia¬ 
tion  sources. 

The  work  described  here  is  a  first  step  in  building  the  WDM  into  a  fully 
configured  radiative-hydrodynamic  description  for  future  load  analysis  and 
design.  We  describe  the  baseline  elements  for  the  method  and  the  successful 
application  of  the  new  tools  to  the  old  and  familiar  problems  with  the  result 
of  some  improved  understanding. 


4.1  Review  of  Fluid  Particle  Methods 

A  fluid  particle  method  may  be  defined  as  any  numerical  technique  that 
replaces  the  conventional  cell  and  boundary  description  of  a  hydrodynamic 
problem  with  some  mobile,  countable  set  of  solution  elements  that  contain 
information  pertinent  to  the  collection  and  not  just  the  mechanical  char¬ 
acteristics  of  an  elementary  fluid  constituent.  Methods  based  on  cellular 
automata  are  then  close  cousins,  being  distinguished  from  fluid  particles  by 
the*  immobile  character  of  the  automata.  So  also  are  finite  element  methods 
on  a  Lagrangian  mesh,  distinguished  by  the  inability  to  mix  and  a  require¬ 
ment  to  specify  the  connectivity  of  the  set.  Fluid  particles  retain  the  grid 
independent  motion  familiar  to  plasma  PIC  methods,  but  add  in  the  physical 
properties  of  the  fluid  they  are  to  represent. 

In  recent  times  fluid  particle  methods  have  been  attempted  by  a  number 
of  workers.  Among  the  most  notable  variations  on  this  theme  we  find. 

Finite  Size  Grid  and  Particle  Hydrodynamics 

The  original  formulation  due  to  Bateson  and  Hewett  [JCP  144,  358- 
378(1998)]  uses  freely  drifting  particles  only  and  focuses  on  particle  annihi¬ 
lation  and  creation  rules  to  track  emerging  features. 

Smooth  Particle  Hydrodynamics 

Developed  originally  by  Monaghan,  it  is  most  common  in  astrophysics 
and  has  been  successfully  specialized  by  Lombardi,  et.al .  [JCP  152,  687- 
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735(1999)]  to  self  gravity  problems.  The  method  uses  PIC-like  particles,  but 
allows  ambient  forces  to  accelerate  them. 

Distribution  Function  Methods 

Generally  not  particle  methods  in  the  strictest  sense,  they  show  a  mix  ol 
Eulerian  and  Lagrangian  technique.  These  methods  can  combine  a  wide  vari¬ 
ety  of  ambient  forces  and  plasma  properties,  e.g.  Huba  and  Lyon  [J.  Plasma 
Physics  61,  391,  1999]  or  Sonendrucker,  et.al.  [JCP  149  687-735(1999)].  The 
distribution  function  is  used  primarily  to  support  flux  and  source  calculations 
for  a  FDTD  advance  of  the  MHD  or  Vlasov  equations. 


4.2  Drift  Kinetic  Grid  and  Particle  Hydrodynamics 

As  a  further  development  of  such  methods  here  we  first  generate  a  second 
order  accurate  solution  of  the  Vlasov  (or  drift  kinetic)  operator  over  a  suffi¬ 
ciently  short  time  interval  that  the  accelerations  and  mean  drifts  are  known 
functions  of  space.  The  new  solution  is  used  as  a  particle  basis,  bringing 
several  advantages  to  the  method  with  respect  to  momentum  transfer  and 
flexibility  in  the  handling  of  external  or  internal  forces. 

Moreover  in  contrast  to  Bateson  and  Hewett,  the  particle  management 
scheme  is  based  upon  a  systematic  energy  and  entropy  conserving  decom¬ 
position  of  observables.  Using  exact  projection  factors  the  exchange  of  in¬ 
formation  from  the  particle  set  to  any  desired  grid  can  be  done  to  machine 
precision.  Such  a  scheme  hence  preserves  any  desired  set  of  observed  initial 
distribution  function  moments  to  forward  projected  moments  of  the  fluid 
particles  that  contain  all  the  information  needed  to  specify  the  source  or  sink 
terms  and  fluxes  for  new  fluid  and  electromagnetic  field  quantities. 

In  contrast  to  earlier  work  we  use  fluid  particles  that  are  a  bit  more 
complex  and  represent  subsets  of  the  fluid  that  may  be  given  any  required 
properties  as  needed,  viz.  temperature,  entropy,  internal  energy,  chemica 
potential,  electrical  conductivity  or  evfen  electrical  charge. 

4.3  Core  Drift  Kinetic  Solution  for  GaPH 

The  primary  strategy  in  the  solution  method  is  an  operator  splitting,  wherein 
the  particles  track  the  free  streaming  processes  and  a  periodic  collection  and 
redecomposition  refines  the  source  and  sink  terms  for  the  external  fields  and 
collisional  processes.  Since  a  physical  system  always  has  both  sets  of  pro¬ 
cesses  going  on  at  once  the  operator  splitting  must  find  a  workable  compro- 
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mise  in  space-time  steps  that  retains  an  image  of  the  system  that  is  “good 
enough”  by  some  set  of  measures. 

The  GaPH  distribution  function  at  any  time  level  is  advanced  from  a 
drifted,  accelerated  Maxwellian  in  a  cartesian  box  h  of  one  to  three  dimen¬ 
sions.  The  minimal  box  dimensions  are  determined  from  the  fundamental 
scale  lengths  for  the  kind  of  fluid  picture  that  is  appropriate  at  any  location 
in  the  problem  domain,  viz.  mean  free  paths,  Debye  lengths,  test  particle 
deflection  distances,  and  the  like.  The  upper  limit  of  a  box  size  is  just  deter¬ 
mined  by  the  desired  resolution  or  the  available  experimental  resolution  for  a 
problem.  The  drift  velocity  C  and  acceleration  A  are  then  Cq°  V?  functions 
of  local  position  x  within  the  box, 

C  =  V  +  x*<5//i,a  =  A+  x*  ajh  . 

Each  box  of  the  problem  is  observed  or  expected  to  contain  a  fixed  initial 
number  of  particles  in  a  known  mix  of  species.  It  is  noteworthy  that  con¬ 
siderable  freedom  exists  in  breaking  up  a  box  of  observables  into  the  fluid 
particles  we  will  use,  because  of  the  presumption  that  the  box  represents  the 
functional  limit  of  resolution  at  any  local  region  of  a  problem. 

The  new  building  blocks  are  by  construction  exact  solutions  of  the 
drift-kinetic  equation  for  the  acceleration  fields  specified,  viz. 


(dt  +  v  •  <9X  +  a  ■  dv)/  =  A)/  =  0  . 


.First  we  can  treat  the  case  of  one  dimensional  particles,  since  multidi¬ 
mensional  ones  are  just  formed  as  outer  products.  To  extract  the  implied 
observables  we  need  the  following  types  of  integrals  over  the  drift  kinetic 
soution  in  the  velocity  subspace. 


N  1 


^h'uVtir 


,  1  ,(((*-  (A  +  ax/h)ty 2  -  x)/t)  -  (V  +  6x/h))^ 
x  exp(— -( - — - )  j 

Integrating  back  in  x  is  moving  back  in  time  along  the  characteristic,  the 
original  Gaussian  velocity  distribution  spreads  in  x  and  v  to  give  hydrody¬ 
namic  motions.  With  accelerations  added,  the  integrals  are  similar  to  those 
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of  Bateson  and  Hewett,  but  now  the  ambient  forces  can  adjust  the  particle 
trajectories  as  they  evolve. 

The  base  particle  solutions  contain  three  common  factors  which  represent 
the  dilation  of  the  box  (due  to  the  shear  in  the  velocity  and  acceleration), 
and  the  (asymmetric)  movement  of  the  box  boundaries: 


.  ,  tS  1  t2a. 

£(f)  =  (1  +  T  +  -ir)  , 


h+(t)  =  h  +  (V  +  <5)  t  +  -  (i4  +  a)  t2  , 


h_(t)  =  -h  +  (V -8)t  +  ^{A-a)t2  . 

The  expected  profile  for  number  density  (per  unit  length,  area,  or  volume) 
n(X,  t)  =<  N  >  1 1  is  then: 


n(X,t)  = 


N 

2hD(t) 


(h+(t)-X)\ 
V2Ut  ) 


+  erf 


((X-h-(t))\\ 

V  y/2  Ut  )) 


while  the  momentum  density  profile  P(X,t)  =<  pv  >  /  is  seen  to  be: 


P^  =  2hDHf)  r  {  ut  )  +  “f  {  V2UT~ ))  (V 


6X  1  ,,  ,  aX 

+  -r  +  2t{A+nr))+ 


mNU  ,(-1/2*^)  _  e(-i/2^#^))  . 
y/2n2hD2(t) 

The  enthalpy  or  total  kinetic  energy  content  profile  H(X,t)  =<  pv 2  >  /2  is 
given  by: 


H(X,t)  = 


mN 

2h£>(t)  X 


(M0z*fl+arf( 

y/2Ut  )  V 


(x-h-(t))\\ 
y/2Ut  )) 


u2  + 


/T,  8X  1  .  .  , 

(y+~r  +  2  *(A+ 

D(t )2 


^))2> 

/ 
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mN_  jj 
+  7^2  hDHty 


V  +  -At  +  X 


(l  6  +  -at^ 


t+2'  h 


J) 


,  ,  (h+(t)-X)2^  ,  (X-h_(t))2, 

(e(-l/2  -W—  )_e(-1/2  — Tg~ )) 


to./V  /i  1  w  /  ,  ,„(fc+<o-x)a,  ,n(x-h_(»))2 

-C/(J  H - (--at)(e^  ^  +  e  ^  , 


V2^2hD*(t)  K  t  '2 


including  both  thermal  and  drift  energies  and  the  changes  in  these  energies 
provided  by  the  external  forces.  There  are  similar  expressions  for  all  velocity 
moments  of  the  core  distribution  function  solution,  such  as  the  heat  flux 
vector  or  the  pressure  tensor. 

Because  the  drift  kinetic  solution  component  is  known  exactly  a  complete 
specification  of  the  influence  that  each  fluid  particle  brings  to  any  other  re¬ 
gion  can  be  tabulated  by  projecting  any  observable  moment  over  any  remote 
domain  of  interest.  The  projection  factors, 


x  exp(  l((((X-(A+ax/Q?/2-T)/t)-(V  +  iX//,))^ 

2  U 

are  -then  the  means  for  collecting  information  out  of  the  particles  to  a  fresh 
set  of  boxes  so  that  the  resolution  of  collisional  processes  and  source  terms 
can  proceed.  In  contrast  to  finite  difference  methods  the  domain  of  support 
can  never  exceed  the  domain  of  influence  for  any  selected  time  step.  More¬ 
over,  summations  of  positive  definite  elements  can  never  produce  a  negative 
particle  or  energy  density. 


4.4  Conserved  Quantities  from  the  Moments 

It  is  obvious  from  the  structures  quoted  above  that  the  moments  of  the  drift 
kinetic  particles  show  conservation  of  number,  momentum,  and  energy  over 
all  space  as  the  particle  spreads  out.  If  we  evaluate  the  various  X>,  t) 

over  the  whole  domain,  the  original  particle  number  is  independent  of  space 
and  time, 
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fh  dx.  ,  1 ,  (((X  —  (A  +  o:x//t)t2/2  x)/t)  (V  +  5x/h)) . 2 


N  r+°°  1  fh  /  dx .  ,  1 ,  ( 

N  =  (2h}  /-oo  dX  UV 2^  ih{T] 6Xp  2  "  U 

The  momentum  ( pV)(X,t )  contains  the  original  momentum  plus  that 
added  by  external  forces, 

1  N  r+°° 

N(5M  +  V)-(j £)/.«» 

1  di.  ,  l,(((X-(A  +  ax//i)t2/2-x)/t)-(F  +  <5x//i))  2 

£VSL(T>  (— )  eXp  ~2  U  ’  '■ 

The  enthalpy  ( pV2/2+pU2/2)(X,t )  contains  the  original  energy  plus  that 
added  by  external  forces  and  force  gradients, 


r). 


1  11  771 TV  /'+00 

(mN/2)  (U2 +  (\t  A  +  V f  +  3(2  ta  +  6 )2)  =  (^)  /  ^  dXx 


/X-x\2  ,  l,(((X-(A  +  ax//i)t2/2-x)A)-(^  +  <5x/h))^ 

/JT>  (— J  eXP('2( - v  >  ' 


4.5  Structure  of  a  Drift  Kinetic  Fluid  Particle 


A  first  example  shows  two  particles  with  an  imbedded  drift  toward  the  right 
and  a  thermal  speed  about  half  this  drift  speed.  Shown  are  the  number 
density  n{X,t),  and  the  enthalpy  density  H(X,t).  Also  shown  are  fields 
derived  from  the  momentum  density  P(X,t).  The  implied  velocity  field, 
V  =  P/mn ,  is  used  to  obtain  the  equivalent  thermal  profile  T(X,t )  = 
H(X,  t)  - 1/2 (p  V2)  and  the  divergence  of  velocity  V  •  V (X,  t).  The  thermal 
profile  has  been  scaled  up  to  illustrate  how  it  blooms  out  relative  to  the 
number  density.  Note  that  the  divergence  of  velocity  spikes  negative  on  the 
leading  side  of  the  particle  where  density  must  increase  as  the  drift  brings 
material  into  a  region.  The  divergence  also  spikes  positive  on  the  trailing 
edge  of  the  particle  as  that  region  rarifies.  The  interior  of  the  particle  retains 
an  essentially  homogeneous  deformation  as  it  drifts.  In  this  manner  the 
continuity  equation  has  been  built  into  the  particle  action  and,  as  shown 
above,  explicitly  conserves  the  total  constituent  number  in  the  drift  kinetic 
fluid  particle  (DKFP). 
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Two  Fluid  Particles  -  the  drift  speed  is  about  twice 


The  “free”  DKFP  accelerates  its  drift  speed  according  to  the  applied 
forces'which  generally  will  include  ambient  drag  from  its  neighbor  particles. 
These  elements  are  free  to  pass  smoothly  through  one  another  as  would 
the  molecules  of  a  gas  in  the  “molecular  flow”  regime.  Note  that  the  time 
origin  and  thermal  speed  are  completely  unconstrained  at  this  level  of  the 
representation.  Each  element  is  a  kind  of  partial  pressure  within  the  gas. 


4.6  Dynamic  Rules  for  the  Momentum 

Not  only  is  number  transport  is  essentially  solved  by  construction,  the  mo¬ 
mentum  transport  is  as  well,  making  it  essentially  redundant  to  ever  compute 
a  pressure  gradient  and  use  that  to  obtain  forces  among  the  DKFP.  The  rea¬ 
son  is  simply  illustrated  in  the  following  figure,  which  portrays  the  velocity 
field  and  the  acceleration  field  (in  addition  to  the  usual  profiles)  while  three 
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Three  Fluid  Particles  Collide  -  the  acceleration  and  velocity  fields  are  everywhere 
proportional,  the  similarity  variable  S(X,t)  is  everywhere  constant. 


DKFP  collide  and  blend  into  one  another.  Each  is  drifting  with  the  same 
(null) ‘velocity  but  subject  to  a  boundary  pressure  gradient  that  drives  them 
together  over  time.  The  similarity  variable  S(X,t )  that  plots  as  a  straight 
line  over  all  space  and  time  is  a  measure  of  the  strictly  proportional*  na¬ 
ture  of  acceleration  to  velocity  during  the  exchange  of  momentum  among  an 
arbitrary  number  “N”  of  DKFP. 

The  fundamental  theorem  concerning  the  momentum  transport  as  de¬ 
picted  above  is  therefore  a  sum,  a  decomposition  of  unity,  viz. 

1  ^  /  tkT _ dxn(^,  t) _ ^ 

~N  -  n(X,t)(V  +  (5X/h)  +  (t/2)(A  +  aX/h))/D(t ))  )i 

showing  that  each  DKFP  expands  and  blends  with  the  sum  of  all  others  in  a 
completely  self-similar  manner  with  the  net  drift  motion  subtracted  out.  So 
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Three  Fluid  Particles  Collide  -  overtime  a  homogeneous  deformation 
results  while  the  similarity  variable  remains  constant. 


right  in  and  requires  no  machinery  to  track  boundaries  or  pressure  gradients 
—  the  boundaries  are  reconstructed  continually  by  the  particle  action. 

The  subsequent  figure  follows  the  same  process  over  time  as  the  colliding 
particles  advance  toward  a  state  of  global  homogeneous  deformation.  Note 
that  the  similarity  variable  remains  strictly  constant  at  one  for  each  particle 
in  the  set  on  all  space  and  time. 

4.7  Isentropic  Rules  for  Decomposition 

In  order  to  represent  a  general  fluid’s  behavior  a  systematic  means  of  dividing 
an  “observables  box"  into  fluid  particles  is  at  least  desireable  and  argueably 
required  if  one  is  to  capture  the  thermodynamics  of  the  system.  Insofar  as 
the  decimation  of  an  “h  box”  into  particles  is  to  be  done  at  a  moment  in 
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time,  as  is  their  recapture  to  any  fresh  set  of  such  boxes  after  a  period  of  drift 
and  interaction,  the  central  thermodynamic  variable  of  interest  is  the  entropy 
of  the  system  defined  by  the  observable  “h  box”.  In  the  figure  is  shown  a  box 
containing  N  =  Nb  +  Nt  base  particles  and  test  particles.  The  test  particles 
are  more  numerous  but  contain  less  material  than  the  base  particles.  The 
spatial  grid  for  the  test  particles  is  distorted  differently  that  for  the  base 
particles.  When  the  box  is  first  divided  into  the  two  classes  at  fixed  volume, 
each  class  is  cooled  to  preserve  the  specific  entropy  which  depends  upon  the 
product  of  number  density  and  temperature  to  the  3/2  power.  Each  class 
can  then  be  distorted  into  spatial  packets  of  variable  size  (smaller  packets 
are  heated  a  bit,  larger  ones  are  cooled)  in  such  a  way  that,  if  the  particles 
were  all  recombined,  the  original  entropy  and  internal  energy  are  recovered 
over  the  whole  box. 

Fluid  Particles  can  decompose  any  observed  set  of  density  or 
velocity  fields  and  preserve  number,  energy  and  entropy. 


(1+ebx) 


(1+eby) 


(1+etx) 


The  fundamental  theorem  concerning  adiabatic  particle  decompositions 
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is  similar  to  that  for  momentum  in  that  it  is  posed  as  a  decomposition  of 
unity.  Summing  over  the  sets  of  base  particles  and  test  particles,  the  various 
distortion  factors  and  etk  are  adjusted  to  preserve  both  the  internal  energy 
and  entropy  of  the  observable  “h  box” .  The  adiabatic  or  isentropic  nature  of 
the  transformation  is  built  in  by  the  power  law  variations  but  the  net  change 
in  internal  energy  (for  an  ideal  gas  in  this  case)  is  a  factor  posed  by  the  left 
hand  side  of 


L,K 

£v«i+' 


,TK.  /_wr_ \ uj 

«w)(l  +%))2/3J  '  \((1  +««)(!  +etj))2/3j 


=  1 


The  decomposition  required  is  achieved  when  a  set  of  distortions  is  found 
such  that  the  equality  is  satisfied.  It  is  easy  to  see  that  if  no  test  particles 
are  selected,  then  the  only  solution  is  for  vanishing  ebk.  Yet  when  some  test 
particles  are  used  the  etk  must  be  chosen  to  boost  their  temperature  back  up 
a  bit  and  balance  the  equation.  The  efficient  multivariate  solution  for  large 
numbers  of  base  and  test  particles  will  clearly  be  a  topic  of  future  effort  but 
problems  of  modest  numbers  are  straightforward.  Note  that  solutions  here 
are  not  unique,  but  link  the  ek  in  a  regular  way.  Note  also  that  the  sum  of 
the  distortions  must  be  zero  since  the  size  of  the  box  remains  invariant  as 
the  base  or  test  particle  parcels  are  compressed  or  stretched. 
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V.  HYDRODYNAMIC  MODELING  OF  NESTED-PLASMA  Z  PINCHES 


I.  Introduction 

Issues  related  to  the  design  and  modeling  of  titanium,  nested-plasma,  Z-pinch  loads  are 
investigated  in  this  report.  Common  to  both  sets  of  issues  is  the  problem  of  determining  initial 
conditions.  They  are  established  in  experiments  nonuniformly  over  the  time  it  takes  for  nested 
loads  to  be  broken  down  and  totally  ionized  as  the  current  rises,  and,  in  computer  calculations  that 
are  used  to  model  these  experiments,  initial  plasma  conditions  are  arbitrarily  chosen  and  taken  to  be 
uniform.  A  second  issue  of  importance  concerns  the  role  of  electron  viscosity  in  diverting  some  of 
the  ion  thermal  energy  into  K-shell  x-rays  in  moderate  atomic  number  Z-pinches.  When  electron 
viscosity  is  ignored  in  modeling  calculations,  all  of  the  viscous  energy  goes  into  ion  heating  which, 
in  turn,  tends  to  divert  more  of  the  compression  heating  into  the  ions  as  well.  However,  the  few 
measurements  of  line  broadening  that  have  been  attempted  to  date  often  tend  to  imply  substantially 
lower  ion  temperatures  than  are  seen  in  calculations. 

Viscous  heating  occurs  when  large  gradients  in  fluid  velocity,  ion  density,  and  electron 
temperature  are  generated  as  shock  fronts  form  during  the  acceleration  and  deceleration  phases 
of  the  pinch  dynamics.  The  resolution  of  such  structures  in  computer  calculations  generally 
requires  many  more  computational  zones  than  present-day  computer  resources  can  handle,  even 
in  one-dimensional  (1-D)  hydrodynamics  calculations.  For  this  reason,  artificial  viscosities  are 
employed  to  model  viscous  heating  since  they  allow  smaller  numbers  of  zones  to  be  used  that 
incompletely  resolve  the  shock  structures  but  still  simulate  viscous  heating  effects.  However,  even 
in  small-zone-number  calculations,  large  gradients  can  be  generated  that  are  only  tens  of  electron 
mean  free  paths  in  length.  In  such  gradients,  other  nonlinear  contributions  to  the  pressure  tensor 
become  important  in  addition  to  the  classical  pressure  tensor  contributions,  which  are  proportional 
to  velocity  gradients.  However,  the  impact  of  these  additional  terms  on  the  pinch  hydrodynamics 
are  not  evaluated  in  this  work. 

The  pinch  modeling  that  is  reported  on  here  utilized  1-D  hydrodynamics  calculations  with  a 
variety  of  initial  conditions  and  zoning  assumptions  in  order  to  determine  some  of  the  sensitivity  of 
the  hydrodynamic  velocity,  temperature,  and  density  gradients  to  these  assumptions.  Since  K-shell 
emissions  depend  on  gradient  structures,  the  dependence  of  these  emissions  on  initial  conditions  is 
ultimately  of  interest.  The  modeling  of  plasma  viscosity  also  presents  a  number  of  problems  since 
(1)  analytic  shock  front  calculations  must  necessarily  ignore  the  ionization  dynamics  and  provide 
limited  benchmarks  against  which  to  check  a  code  calculation,  (2)  real  viscosities  are  not  used 
in  computer  codes  for  lack  of  zoning  resolution,  (3)  electron  viscosities  only  become  important 
beyond  certain  levels  of  ionization  in  moderate  to  high  atomic  number  plasma  and  then  only 
under  certain  plasma  conditions,  and  (4)  accurate  computer  modeling  of  shock  front  generation  in 
moderate  atomic  number  plasmas  requires  good  ionization  dynamic  and  radiation  transport  models 
in  addition  to  a  large  number  of  zones.  Consequently,  for  the  modeling  described  in  this  report, 
the  standard  form  of  artificial  viscosity  was  used,  but  the  heating  it  produced  was  divided  between 
electrons  and  ions  according  to  the  ratio  of  their  classical  viscosities,  a  ratio  that  was  recently 
computed  analytically  [K.  G.  Whitney,  “Momentum  and  Heat  Conduction  in  Highly  Ionizable 
Plasmas”,  Phys.  of  Plasmas,  6,  816  (1999)]. 
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Two  sets  of  calculations  of  titanium  and  krypton  plasma  implosions  were  carried  out.  In 
one  set  of  calculations,  the  imploded  loads  were  initially  uniform  fills  of  krypton  and  titanium 
plasmas  at  a  temperature  of  4  eV.  These  calculations  were  carried  out  using  60  zones,  and  the  fills 
were  driven  by  a  linearly  rising  current  pulse  that  was  terminated  when  a  prescibed  energy  input 
was  achieved.  The  mass  of  the  fill,  its  initial  density,  and  the  energy  input  were  varied  so  that 
comparisons  of  the  ion  temperature  profiles  and  K-shell  yields  and  powers  could  be  made  with  and 
without  electron  viscosity  in  the  calculations.  Fill  masses  between  10  and  1.25  mg/cm  were  used. 
A  collisional-radiative  equilibrium  (CRE)  ionization  dynamics  was  employed;  ultimately  however, 
more  accurate  modeling  of  shock  fronts  will  require  time-dependent  ionization  since  the  passage 
of  a  shock  causes  rapid  heating  to  occur  that  drives  the  ionization  out  of  equilibrium.  The  results  of 
the  uniform-fill  calculations  are  described  in  section  II.  Next,  a  nested  titanium  plasma  geometry 
was  set  up  in  which  an  an  outer  portion  of  the  above  described  fill  plasma  was  separated  from  the 
inner  portion  by  moving  the  outer  plasma  out  in  radius,  thereby  creating  a  small  mass  of  reduced 
density  (by  a  factor  of  10)  between  the  inner  and  outer  plasmas.  This  nested  plasma  load  was 
again  driven  by  a  linearly  rising  current  to  pre-determined  energy  inputs.  A  mass  of  1.25  mg/cm 
was  used  in  a  number  of  nested-plasma  calculations,  since  it  corresponds  to  an  amount  of  plasma 
mass  (and  not  necessarily  wire  mass)  that  can  be  imploded  on  the  Z  accelerator  to  velocities  that 
are  needed  to  effectively  excite  the  K-shell  of  titanium.  The  results  of  these  calculations,  with  and 
without  electron  viscosity,  are  described  in  section  III.  Finally,  in  section  IV,  the  K-shell  power  and 
yield  behavior  of  a  series  of  calculations  that  correspond  in  essential  ways  to  a  series  of  titanium 
nested  array  shots  that  were  conducted  on  the  Z  accelerator  are  discussed.  These  calculations  use 
the  same  geometry  as  the  section  HI  nested  plasma  loads.  The  experiments  were  initiated  using 
two  sets  of  nested  titanium  wire  arrays  of  varying  total  mass  with  a  2:1  mass  ratio  and  a  2:1  radii 
ratio  of  outer  to  inner  plasmas.  These  loads  were  accelerated  to  larger  implosion  velocities  as 
the  load  mass  was  reduced.  The  nested  plasma  calculations  described  in  section  IV  were  carried 
out  using  the  same  masses  and  achieving  approximately  the  same  implosion  velocities  as  in  the 
Z  experiments.  Comparisons  between  these  experiments  and  the  corresponding  calculations  and 
some  comments  about  future  work  are  made  in  section  V. 

II.  Analytic  Considerations 

When  B-field  effects  are  ignored,  the  classical  electron  and  ion  pressure  tensors,  ne  and  Ir¬ 
respectively,  are  linearly  proportional  to  the  symmetrized  and  traceless  version,  (Vv/)s,  of  the 
velocity  gradient  tensor: 


De  =  i?eo(Vv/)s,  Hi  =  r/io(Vv/)s. 


Thus,  ne/Hi  =  r  =  T)eQ/7iio,  and,  from  the  calculation  of  7]e0  and  /rji o, 

1.33Z(Z2  +  2.820Z  +  1.343)  /  me  \  1/2  /Te  \  5/2 

r~Z3  +  4.434 Z2  +  5.534 Z  +  1.780  \Mion)  \Tt ) 

This  ratio  is  a  function  only  of  the  degree  of  plasma  ionization,  Z,  and  the  electron  to  ion  temperature 
ratio.  Also,  since 

ne  -|-  Hi  =  Harit/tcta/ j 
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one  can  partition  the  artificial  viscous  heating  between  the  electrons  and  ions  using  the  above  easily 
calculable  ratio,  r.  This  procedure  represents  the  minimum  step  that  can  be  taken  in  evaluating  the 

role  of  electron  viscosity  in  Z-pinch  implosions. 

The  role  expected  can  be  discerned  from  an  examination  of  the  fluid  equations  for  the 
electron  and  ion  temperatures.  These  equations  have  similar  forms  when  written  in  terms  of 
the  dimensionless  viscosity  coefficients,  rjeo  and  rji0,  which  are  defined  by 

Veo  =  rjeO  PeTe,  TjiQ  =  ViO  PiTi- 

When  written  in  terms  of  these  coefficients,  the  temperature  equations  become 

(a,  +  T.  +  |  (~  («>,)  -  for.  (vvy)  T,  =e-i  heating  & 

cooling  terms, 

(a,  +  vX)Ti  +  l (;|;(r«r)  -  w(vv/)!)s  =  e-i cooling, 

Both  congressional  and  viscous  heating  enter  these  equations  as  variable  exponential  growth 
rates.  Additional  heating  and  cooling  takes  place  locally  due  to  electron-ion  collisions  and  thermal 
conduction,  while  global  cooling  occurs  because  of  radiation  losses.  These  equations  show  that 
plasma  compression  acts  to  heat  both  electrons  and  ions  at  the  same  rate  in  the  absence  of  viscous 
heating.  However,  when  electron  viscosity  is  ignored  and  all  of  the  viscous  heating  is  put  into  the 
ions,  the  exponential  rate  of  ion  heating  exceeds  the  electron  rate,  which,  in  this  case,  is  due  to 
compression  heating  alone.  Thus,  the  ions  become  significantly  hotter  than  the  electrons.  This  effect 
shows  up  in  all  computer  code  calculations  that  include  ion  viscosity  only.  The  addition  of  electron 
viscosity  should,  therefore,  reduce  ion  temperatures  and  raise  electron  temperatures  slightly,  which 
suggests  the  possibility  of  converting  some  of  this  energy  into  extra  K-shell  emission.  However, 
this  possibility  will  depend  on  the  plasma  conditions  that  are  generated  during  the  implosion  since 
the  ratio,  r,  is  substantially  reduced  when  T»  is  allowed  to  become  much  larger  than  Te.  Radiation 
cooling  can  become  a  drag  on  the  electron  temperature  that  acts  on  the  ions  only  indirectly  through 
electron-ion  collisions.  Thus,  the  ions  can  get  hotter  than  the  electrons  when  radiation  losses  are 
high  and  the  electron-ion  coupling  is  too  weak.  Under  conditions  of  too  low  an  ion  density,  the 
rerouting  of  ion  thermal  energy  into  K-shell  emissions  is  inhibited. 

III.  Uniform  Fill  Calculations 

Since  the  ratio  of  electron  to  ion  viscosity  increases  roughly  as  Zs/2  for  K-shell  emitters, 
two  uniform-fill  krypton  calculations  were  carried  out  prior  to  uniform-fill  titanium  calculations  to 
test  for  the  influence  of  this  Z-scaling.  The  mass  of  the  krypton  load  was  taken  to  be  10  mg/cm, 
its  initial  radius  was  4  cm,  and  the  energy  input  to  the  pinch  corresponded  approximately  to  an 
77*  of  two.  Results  from  these  calculations  are  shown  in  Figs.  (l)-(5).  Because  60  zones  were 
used  in  the  calculations,  the  artificial  viscosity  was  not  enhanced.  However,  even  with  only  60 
zones,  very  sharp  gradients  in  temperature,  density,  fluid  velocity,  and  fractional  viscous  heating 
were  generated,  as  illustrated  in  Figs.  (1)  and  (2)  at  two  different  times  (60  and  80  ns)  early  in 
the  implosion.  Figs.  (3)-(5)  show  comparisons  between  CRE  calculations  that  were  done  wi 
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and  without  electron  viscosity.  When  included  in  the  calculation,  electron  viscosity  produced 
significant  reductions  in  the  ion  temperature:  illustrated  in  the  ion  temperature  profiles  of  Figs. 
(3)  and  (4)  at  80  and  113.5  ns  respectively.  In  this  particular  calculation,  reduced  ion  heating  led 
to  a  20%  increase  in  K-shell  energy  output  and  produced  close  to  4  MJ.  As  Fig.  (4)  shows,  all 
of  the  K-shell  emission  originated  from  the  core  of  the  compressed  plasma.  The  outer  portion  of 
plasma,  which  contains  a  sizable  amount  of  the  plasma  mass,  emitted  only  in  the  L-shell.  Prior  to 
K-shell  emission,  the  electron  temperature  profile  (at  80  ns  in  Fig.  (3))  shows  only  slight  changes 
due  to  electron  viscosity  in  contrast  to  the  ion  temperature  profile  since  there  are  Z  times  as  many 
electrons  to  be  heated  as  ions. 

When  10  mg/cm  of  a  uniform- fill  of  titanium  (at  an  initial  radius  of  2  cm)  was  substituted 
for  the  10  mg/cm  of  uniform-fill  krypton,  similar,  but  smaller,  reductions  in  ion  temperature  were 
seen.  These  reductions  in  the  ion  thermal  energy  again  showed  up  as  increases  in  K-shell  emission. 
Figs.  (6)-(ll)  illustrate  this  behavior.  Figs.(6)-(8)  contain  histories  of  the  K-shell  power  and 
energy  outputs  in  three  calculations  carried  out  with  and  without  electron  viscosity.  Figs.  (9)-(ll) 
contain  profiles  of  the  ion  temperature  at  times  just  before  peak  compression  (with  and  without 
electron  viscosity)  and  of  the  average  degree  of  ionization.  The  energy  input  to  the  plasma  was 
systematically  increased  in  these  calculations.  In  Figs.  (6)  and  (9),  the  energy  input  corresponded 
to  r]*  =  1.8,  in  Figs.  (7)  and  (10),  this  input  was  doubled  to  77*  =  3.6,  and,  in  Figs.(8)  and  (11),  it 
was  doubled  again  to  77*  =  7.2.  These  doublings  were  accomplished  by  reducing  the  risetime  of 
the  current. 

The  results  in  Figs,  (l)-(l  1)  were  not  obtained  through  a  systematic  and  thorough  study  of  how 
electron  viscosity  affects  moderate  atomic  number  pinch  implosions;  hence,  they  only  represent  a 
first  look  at  how  the  electron  viscosity  can  influence  K-shell  emission  and  about  how  this  influence 
depends  on  plasma  conditions  and  on  the  atomic  number  of  the  plasma.  Figs.  (6)-(8)  demonstrate 
that  the  percentage  increase  in  K-shell  output  due  to  reductions  in  ion  temperature  will  vary  with 
plasma  conditions  and  that,  for  the  fill  plasma  chosen,  it  exhibits  a  maximum  as  a  function  of 
the  input  energy.  Figs.  (5)-(8)  demonstrate  that  this  percentage  can  go  up  as  a  function  of  atomic 
number  as  the  Z-scaling  of  the  ratio  of  viscosity  coefficients  suggests.  Note  also  that,  as  the  energy 
input  in  these  titanium  calulations  was  doubled,  the  core  ion  temperatures  more  than  doubled,  the 
core  ionization  levels  approached  K-shell  bum-through,  and,  as  a  consequence,  the  K-shell  yield 
at  first  doubled  but  then  did  not. 

Since  the  Z-accelerator  does  not  have  the  energy  to  drive  10  mg/cm  loads,  two  sets  of  uniform- 
fill  calculations  were  carried  out  to  observe  what  happens  when  the  imploded  mass  is  progressively 
lowered.  The  mass  of  the  titanium  load  was  lowered  by  first  increasing  the  initial  density  of  the 
load  by  \/2  and  then  by  reducing  the  initial  radius  of  the  load,  preparatory  to  doing  nested-plasma 
calculations.  In  one  set,  the  load  mass  was  reduced  to  5  mg/cm  (for  77*  =  3.6)  and,  in  the  other,  to 
a  mass  of  1 .25  mg/cm  for  an  77*  of  2,  which  corresponds  to  a  mass  and  energy  input  obtainable  in 
Z-accelerator  experiments.  Results  from  these  calculations  are  shown  in  Figs.  (12)-(15).  In  each 
set  of  calculations,  the  inclusion  of  electron  viscosity  reduced  the  ion  temperature  [Figs.  (14)  and 
(15)].  Note  that,  by  halfing  the  mass  of  the  load  and,  thus,  the  number  of  radiators,  the  K-shell 
energy  output  was  halfed  even  though  the  power  output  was  not  significantly  reduced  [compare 
Fig.  (12)  with  Fig.  (7)].  The  energy  radiated  in  the  1.25  mg/cm  calculation  [Fig.  (13)]  also  went 
down  in  proportion  to  both  the  number  of  radiators  and  the  input  energy.  Figs.  (14)  and  (15)  indicate 
a  tendency  for  the  plasma  core  to  be  over-ionized  in  these  calculations,  i.e.,  to  be  non-optimally 
compressed.  Note,  both  in  Figs.  (6)-(ll)  and  in  Figs.  (12)-(15),  that  ion  temperatures  generated 
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in  the  core  region  of  the  plasma  increased  considerably  as  the  amount  of  K-shell  emission  that  the 
plasma  produced  from  this  region  increased.  Thus,  even  though  the  electron  viscosity  brought  ion 
temperatures  down  in  these  particular  calculations,  it  could  not  keep  them  down.  This  tendency 
for  the  ions  to  overheat  is  detrimental  to  K-shell  emission  since  it  reduces  plasma  compression  and 
radiation  rates,  and  it  stores  energy  in  the  ions  that  cannot  be  radiated  away. 


IV.  Nested-plasma  Calculations 

The  1.25  mg/cm  uniform-fill  calculation  was  used  as  a  stepping  stone  to  the  1.25  mg/cm 
nested-plasma  calculations  that  are  described  in  this  section.  They,  in  turn,  provided  guidance  for 
the  five  nested-plasma  calculations  that  were  made  in  correspondence  with  five  nested-titanium, 
Z-accelerator,  wire  shots  and  that  are  described  in  the  following  section.  Several  of  the  1 .25  mg/cm 
calculations  were  made  to  evaluate  the  way  zoning,  mass  distribution,  and  implosion  velocity 
affect  the  formation  of  plasma  gradients  and  their  influence  on  x-ray  outputs.  The  influence  of 
electron  viscosity  on  nested  titanium-plasma  implosions  was  also  investigated.  The  nesting  used 
in  the  calculations  was  done  by  dividing  the  1.25  mg/cm  uniform-fill  plasma  into  two  parts  and 
by  moving  the  outer  part  out  in  radius,  by  roughly  a  factor  of  two.  Since  one  cannot  have  a 
perfect  vacuum  between  plasmas  in  a  computer  calculation,  the  calculations  were  zoned  so  that  the 
density  between  the  inner  and  outer  plasmas  could  be  reduced  by  a  factor  of  ten.  Thus,  the  nesting 
used  in  these  calculations  only  approximated  certain  of  the  essential  features  of  the  Z-accelerator 
experiments,  and  for  the  following  additional  reasons.  One,  the  calculations  described  here  are 
of  the  shell-on-fill  rather  than  shell-on-shell  variety.  To  remove  plasma  from  the  axis  for  shell- 
on-shell  implosions  requires,  in  principle,  even  greater  numbers  of  zones  and  more  computer  time 
than  were  used  in  the  shell-on-fill  calculations.  Two,  the  initial  plasma  conditions  that  must  be 
chosen  when  modeling  the  experiments  are  essentially  unknown  and,  at  best,  approximate.  Shell- 
on-shell  calculations  have  an  addition  parameter  to  play  with,  since  the  initial  shell  thicknesses 
are  chosen  somewhat  arbitrarily.  Three,  the  plasma  densities  generated  early  in  the  experiments 
are  not  only  unknown  but  also  non-uniform.  In  addition,  the  state  of  the  inner  titanium  (wires 
or  plasma)  at  the  time  of  impact  is  unknown.  Altogether,  these  problems  imply  that  there  is  no 
rational  for  attempting  calculations  that  purportively  accurately  simulate  (rather  than  correspond) 
to  current  nested-wire  experiments.  The  shell-on-fill  calculations  described  here  have  the  utility 
of  demonstrating  hydrodynamic  behavior  that  promotes  K-shell  emission  from  moderate  atomic 
number  loads.  Another  of  their  benefits  is  that  they  suggest  certain  hydrodynamic  behavior  that 
might  favor  wire-on-foam experiments  with  pre-ionization  over  simple  wire-on-wire  experiments. 

The  seven  cases  listed  in  Table  I  were  analyzed  and  compared  in  these  initial  calculations. 

In  the  first  three  cases,  the  initial  plasma  was  kept  unchanged,  but  the  implosion  was  calculated 
employing  different  zoning  schemes.  In  case  1,  only  80  zones  were  used,  16  of  which  were 
allocated  to  the  inner  (nested)  plasma.  In  cases  2  and  3,  the  total  number  of  zones  was  increased 
to  100,  and  the  number  used  to  resolve  the  gradients  generated  in  the  inner  plasma  was  17  m  case 
2  and  27  in  case  3.  The  dynamics  of  these  nested-plasma  implosions  is  illustrated  in  Fig.  (16).  it 
contains  time  histories  of  the  shell  boundaries  of  the  outer  plasma  (in  red)  as  they  implode  inward 
and  impact  the  outer  boundary  of  the  nested  plasma  (in  blue)  at  approximately  46.5  ns.  The  second 
blue  curve  in  Fig.  (16)  shows  the  time  history  of  the  first  cell  boundary.  Even  though  a  perfect 
vacuum  does  not  exist  between  the  two  plasmas  in  the  calculation,  the  collision  time  is  still  fairly 
well  defined  Note  that  the  outer  surface  of  the  outer  plasma  bounces  outward  at  the  time  when  the 
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Table  I.  Nested-plasma  Calculations 


Case 

# 

#  Zones 
(Total) 

#  Inner 
Plasma  Zones 

Mass  Ratio 
(Outer/Inner) 

Inner  Plasma 
Radius  (mm) 

V* 

1- 

80 

16 

2:1 

4.2 

2 

2. 

100 

17 

2:1 

4.2 

2 

3. 

100 

27 

2:1 

4.2 

2 

4. 

100 

35 

1:1 

4.3 

2 

5. 

100 

27 

2:1 

4.2 

4 

6. 

100 

31 

2:1 

2.2 

2 

7. 

100 

35 

1:1 

4.3 

4 

current  drive  is  terminated.  Three  vertical  lines  are  drawn  in  Fig.  (16)  at  47.8,  49.3,  and  50.7  ns. 
They  identify  the  three  times  in  the  implosion  history  at  which  the  spatial  profiles  of  temperature, 
density,  and  ionization  are  drawn  in  Figs.  (17)  and  (18).  These  two  figures  illustrate  that  strong 
shocks  are  formed  during  the  collision  of  the  outer  and  inner  plasmas.  However,  even  within  a 
100  zone  calculation,  these  shocks  are  poorly  resolved.  The  profiles  in  Fig.  (17)  show  how  the 
temperature  front  moves  ahead  of  the  density  front,  leading,  at  the  time  of  peak  compression,  to 
a  less  dense  core  plasma.  Since  K-shell  x  rays  originate  from  this  core  plasma,  any  reductions  in 
its  density  can  seriously  impact  this  emission.  Fig.  (18)  shows  that  during  shock  formation  the 
plasma  is  driven  into  the  ground  helium-like  state  of  the  K-shell,  but  that  [Fig.  (19)]  actual  K-shell 
emission  does  not  begin  until  52  ns  just  before  the  final  compression  and  heating  take  place. 

Two  case  3  calculations  were  made,  one  with  and  one  without  electron  viscosity.  The  addition 
of  electron  viscosity  allowed  an  additional  20  kJ  of  K-shell  x  rays  to  be  extracted  from  the  plasma 
[Fig.  (19)].  Fig.  (20)  shows  the  reason  for  this  added  emission.  Two  ion  temperature  profiles  are 
drawn  at  the  onset  of  K-shell  emission.  Without  electron  viscosity,  the  ions  are  much  hotter  than 
the  electrons  in  the  core  region  of  the  plasma  as  the  emission  begins.  Consequently,  the  temperature 
ratio  in  the  formula  for  T]e0 /r}i0  is  reduced  allowing  some  of  the  final  ion  heating  to  be  averted  and 
channeled  instead  into  x  rays. 

Another  variable  that  can  be  used  in  place  of  the  radius  when  plotting  profiles  is  m(r),  the 
plasma  mass  contained  within  a  cylinder  of  radius,  r : 

m(r)  =  2x  /  dr'r' p(r'). 

Jo 

When  <  Z  >  is  plotted  as  a  function  of  m(r)  rather  than  of  r  at  the  time,  t  —  51.6,  just  as  K-shell 
emission  begins,  one  learns  that  only  ~  28%  of  the  plasma  mass  in  this  load  has  been  primed  to 
emit  from  the  K-shell  as  the  final  compression  stage  of  the  implosion  begins. 

Cases  1-3  address  the  question,  how  does  the  number  and  location  of  the  hydrodynamic  zones 
in  a  computer  calculation  affect  the  plasma  dynamics  and  plasma  profiles  that  are  calculated?  Part 
of  the  answer  is  given  in  Figs.  (21)-(23).  The  profiles  shown  in  these  figures  were  all  calculated  at 
the  time  of  peak  compression.  The  plasma  compression  was  significantly  reduced  along  with  the 
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core  density  by  going,  in  case  1  to  case  2,  from  80  to  100  zones,  while  the  core  electron  temperature 
was  increased,  and  the  ratio  of  electron-to-ion  temperature  was  reduced.  The  increased  zoning 
of  the  inner  plasma,  in  going  from  case  2  to  case  3,  led  to  increased  gradient  resolution  in  die 
core,  which  increased  the  core  compression  without  changing  its  temperature  or  its  electron-to-ion 
temperature  ratio.  Figs.  (21)-(23)  provide  examples  of  how  the  computed  hydrodynamics,  which  is 
driven  by  gradients,  is  a  function  of  the  zoning  used  for  small  zone  number,  since  zoning  determines 

the  gradient  structure  of  the  plasma  that  can  be  computed. 

Case  4  addressed  the  question,  could  mass  shifting  and  larger,  initial,  inner-plasma  densities 
increase  the  mass  participation  from  the  ~  28%  level  seen  in  Fig.  (20)?  In  case  4,  the  increased 
inner  density  was  achieved  by  shifting  mass  from  the  outer  to  the  inner  plasma  and  making  the 
masses  equal.  In  case  4,  the  plasma  received  the  same  input  energy  as  case  3.  Figs.  (24)  and  (25) 
show  that  this  mass  and  density  shift  did  increase  the  amount  of  K-shell  mass  participation  (drawn 
in  Fig.  (24)  at  the  times  of  peak  compression)  as  well  as  the  K-shell  yield  by  increasing  the  peak 
power  output.  Notice,  however,  that  the  power  pulse  was  clipped  following  peak  power.  Thus,  the 
energy  output  was  limited  for  this  particular  load  because,  at  tj  =2,  there  was  no  more  energy 
available  to  be  radiated  from  the  core,  i.e.,  the  mass  participation  and  energy  content  of  the  core 
in  this  case  did  not  increase  enough.  Both  Figs.  (19)  and  (25)  contain  a  plot  of  the  outer  radius 
of  the  inner  plasma.  These  curves  convey  information  about  the  elasticity  or  inelasticity  of  the 
compressions  and  suggest  that  the  L-shell  radiation  losses  in  case  4  were  larger  than  those  in  case 
3  (as  they  were  by  approximately  100  kJ). 

Case  5  looked  at  the  response  of  the  case  3  plasma  when  the  current  was  increased  and  the 
implosion  was  driven  by  twice  the  input  energy  of  case  3.  For  this  plasma  configuration,  however, 
a  doubling  of  the  input  energy  did  not  double  the  K-shell  power  or  the  energy  output  [see  Fig.  (26)]. 
The  L-shell  output  for  this  case,  however,  was  less  than  that  of  case  4.  Note  that  the  energy  output 
in  these  calculations  changes  in  proportion  to  the  peak  power  output  generated  at  the  time  of  peak 
compression  since  the  current  is  turned  off  prior  to  peak  compression.  At  the  start  of  the  x-ray 
pulse  in  case  5,  when  25%  of  the  plasma  is  already  beginning  to  radiate  in  the  K-shell,  the  electron 
viscosity  had  acted  to  reduce  the  ion  temperature  [Fig.  (27)];  however,  the  subsequently  large 
x-ray  output  keeps  the  electrons  cool  as  the  ions  are  heated,  resulting  in  little  difference  between 
the  K-shell  outputs  with  and  without  electron  viscosity  [Fig.  (26)].  The  increased  energy  output 
occurred  as  a  result  of  the  approximately  14%  increase  in  mass  participation  that  the  extra  energy 
input  was  able  to  promote.  Fig.  (28)  demonstrates  this  result  by  comparing  ionization  levels  as 
a  function  of  m(r)  at  the  time  of  peak  emission  for  the  case  3  and  5  plasmas.  Fig.  (29),  when 
compared  to  Fig.  (18),  explains  why  this  increase  occurred.  The  added  energy  input  puts  more  of 
the  plasma  into  the  helium-like  ground  state  early  on  in  the  implosion,  preparatory  to  its  lighting-up 
during  the  final  stages  of  compression. 

The  last  two  cases  studied,  6  and  7,  are  illustrated  in  Figs.  (30)-(32).  In  case  6  in  contrast  to 
case  4,  the  mass  ratio  of  outer-to-inner  plasma  was  kept  at  2: 1,  while  the  initial  density  of  the  inner 
plasma  was  raised  by  decreasing  its  radius.  The  effect  of  increasing  density  by  decreasing  radius 
was  to  increase  the  power  output  (to  that  of  case  5)  but  to  decrease  the  energy  output  by  clipping 
the  power  pulse  even  more  severely  than  in  case  4.  That  this  clipping  effect  could  be  mitigated 
was  demonstrated  in  case  7.  The  plasma  in  case  7  was  the  same  as  in  case  4,  but  the  energy  input 
was  doubled.  Hence,  case  7  is  to  case  4  as  case  5  is  to  case  3.  Figs.  (19)  and  (26)  show  that 
doubling  the  energy  input  when  the  outer-to-inner  plasma  mass  ratio  is  2:1  does  not  double  the 
energy  output,  while  Figs.  (25)  and  (31)  show  that  doubling  the  energy  input  more  than  doubles 
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the  output  when  the  mass  ratio  is  1:1.  Again,  the  large  electron  cooling  that  the  K-shell  output 
produced  eliminated  the  possibility,  in  this  case,  of  seeing  extra  K-shell  output  as  a  consequence  of 
electron  viscosity  [Fig.  (31)],  even  though  a  large  difference  in  ion  temperatures  at  the  beginning  of 
the  K-shell  pulse  existed  in  the  two  calculations  (with  and  without  electron  viscosity)  [Fig.  (32)]. 
Note  also  in  Fig.  (32)  that  the  doubled  energy  input  in  case  7  increased  K-shell  mass  participation 
to  ~  50%. 

V.  Variable  Mass,  Variable  77*,  Nested  Calculations 

Five  different  nested-titanium-wire  experiments  were  carried  out  on  the  Z  accelerator.  For 
two  of  the  five,  3  shots  each  were  done.  Table  II  contains  the  basic  parameters  of  the  9  shots. 

Table  II.  Experimental  Parameters 


Shot 

# 

Outer  Diam 
(mm) 

Outer  Mass 
Hg/cm 

Inner  Mass 
(Mg/cm 

Wire  Gap 
Spacing  (mm) 

V 

422 

40 

1246 

623 

0.9 

1.6 

BBS! 

45 

534 

1.18 

2.1 

45 

1072 

534 

1.18 

2.1 

468 

45 

1072 

534 

1.18 

2.1 

424 

858 

429 

1.64 

2.6 

465 

50 

858 

429 

1.64 

2.6 

471 

858 

429 

1.64 

2.6 

470 

55 

714 

357 

2.16 

3.1 

467 

60 

500 

250 

3.37 

3.8 

As  mentioned  earlier,  two  problems  immediately  arise  in  attempting  to  do  calculations  that  closely 
simulate  these  experiments.  One,  different  wire  array  diameters  were  used  in  the  experiments, 
which  produced  large  variations  in  wire  gap  spacing:  from  0.9  mm  in  shot  422  to  3.37  in  shot  467. 
The  effect  of  gap  spacings  on  the  implosion  dynamics  cannot  be  modeled,  but  their  effect  in  the 
experiments  is  known  to  be  large.  This  problem  arises  because  wire  number  cannot  be  increased 
indefinitely  on  a  given  machine  as  the  array  diameter  is  increased.  An  analogous  theoretical  zoning 
problem  arises  because  the  number  of  zones  in  a  computer  calculation  must  also  be  increased  as 
the  load  radius  is  increased;  however,  the  number  of  zones  cannot  be  increased  indefinitely  due 
to  limitations  in  computer  resources.  Two,  because  the  wire  explosion  and  merging  dynamics  is 
unknown,  the  initial  conditions  to  be  used  in  the  computers  calculations  are  basically  unknown.  In 
either  1-D  or  2-D  calculations,  one  can  only  roughly  approximate  them.  For  these  two  reasons, 
the  calculations  that  were  done  were  chosen  to  correspond  with  the  experiments  in  the  following 
way.  Plasma-on-plasma  collisions  were  modeled  in  a  fixed  geometry  with  roughly  a  2: 1  ratio  in 
outer-to-inner  plasma  diameters  in  accord  with  the  experiments.  Also,  the  same  masses  were  used 
as  in  the  experiments,  and  they  were  imploded  to  roughly  the  same  velocities  as  in  the  experiments. 
Results  from  these  calculations  are  shown  in  Table  m  and  in  Figs.  (33)-(47). 


106 


Table  III.  Experimental/Calculational  Comparisons 


Shot 

# 

K-shell 

Yield 

(kJ) 

Total 

Yield 

(kJ) 

K-Mass 

Radiating 

(mg/cm) 

Mass 

Fract 

% 

Exp  K-shell 
Yield,  Main 
Pulse  (kJ) 

Exp  Total 
Yield,  Main 
Pulse  (kJ) 

422 

188 

1005 

0.35 

18.7 

47.2 

700 

423 

254 

1037 

0.37 

23.1 

54.5 

672 

466 

254 

1037 

0.37 

23.1 

41 

620 

468 

254 

1037 

0.37 

23.1 

45 

598 

ave 

254 

1037 

0.37 

23.1 

46.8 

630 

424 

269 

870 

0.36 

27.9 

56.2 

578 

465 

269 

870 

0.36 

27.9 

55 

874 

471 

269 

870 

0.36 

27.9 

68 

650 

ave 

269 

870 

0.36 

27.9 

59.7 

701 

470 

249 

708 

0.33 

30.8 

84 

626 

467 

140 

269 

0.27 

35.3 

57 

364 

The  first  five  figures,  Figs.  (33)-(37),  show  calculated  K-shell  power  and  energy  outputs  as  a 
function  of  time  as  the  load  mass  was  progressively  decreased  and  the  load  input  energy-per-ion  and 
implosion  velocity  was  increased.  The  powers  and  yields  at  first  increase  achieving  a  maximum  for 
the  1.29  mg/cm  (77*  =  2.6)  load  before  falling  again.  This  behavior  is  explained  in  Figs.  (38)-(42). 
These  figures  contain  plots  of  the  ionization  level  of  the  plasma  at  the  time  of  peak  emission  as 
a  function  of  the  mass  enclosed  in  cylinders  of  radii,  r.  In  all  cases,  the  K-shell  emission  comes 
from  the  core  of  the  plasma,  which  grows,  at  first,  slightly  in  size  and  then  diminishes  slightly.  The 
temperature  and  density  profiles  at  the  time  of  peak  compression  and  emission  (in  Figs.  (43)-(47)) 
explain  the  behavior  seen  in  Figs.  (38)-(42).  As  more  energy-per-ion  is  imparted  to  the  load,  core 
temperatures  are  driven  higher,  leading  to  core  overheating.  Thus,  core  densities  go  down,  reducing 
x-ray  radiation  rates.  This  behavior  is  indicative  of  strong  shock  wave  convergence  on  axis. 

VI.  Conclusions 

One  of  the  conclusions  that  can  tentatively  be  drawn  from  the  above  comparisons  is  that 
substantial  improvements  in  load  performance  are  possible  through  load  designs  that  better 
approximate  the  initial  conditions  used  in  the  calculations.  Of  course,  the  calculations  assume 
overly  ideal  initial  conditions,  i.e.,  they  assume  that  plasmas  are  preformed  and  uniform  in  both 
density  and  temperature.  Consequently,  the  calculations  produce  much  larger  K-shell  energy 
outputs  in  a  single  high-power  pulse  than  the  experiments  do.  Experiments,  on  the  other  hand, 
are  carried  out  in  ways  that  make  it  impossible  to  create  such  ideal  plasma  initial  conditions. 
Consequently,  it  is  not  entirely  surprising  that  they  produce  up  to  a  factor  of  five  smaller  K-shell 
outputs  in  the  main  pulse  than  calculations  do.  Some  of  this  large  difference  can  be  attributed  to  the 
poorly  understood  and  uncontrolled  explosion  dynamics  of  the  wires.  However,  the  calculations 
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suggest  that  generating  more  uniform  plasmas  initially  could  produce  large  payoffs  in  significantly 
improved  K-shell  x-ray  outputs  from  moderate  atomic  number  Z-pinches.  Under  the  assumption 
that  approximate  agreement  between  trends  seen  in  the  theoretical  and  experimental  implosions 
is  due  to  common  features  in  their  dynamics  (even  when  their  initial  conditions  have  marked 
differences),  one  can  also  make  tentative  interpretations  of  the  experimental  data  and,  consequently, 
some  proposals  for  experiments  that  could  improve  load  performance. 

One  notable  trend  common  to  both  the  experiments  and  calculations  is  the  appearance  of  a 
K-shell  yield  maximum.  The  calculations  show  that,  as  the  load  mass  is  lowered  and  the  input 
energy  or  implosion  velocity  is  raised,  the  core  region  is  driven  to  higher  temperatures  and  lower 
densities.  In  collisional-radiative  equilibrium,  which  is  assumed  in  the  calculations,  these  trends 
lead  to  K-shell  bum-through  and  lower  power  emissions,  which  may  or  may  not  be  present  in  the 
experiments.  In  addition,  as  the  load  mass  was  reduced  in  the  calculations,  the  amount  of  mass 
in  the  core  that  ended  up  radiating  in  the  K-shell  went  slowly  downward  at  first  and  then  more 
rapidly  downward  (even  as  the  mass  fraction  went  steadily  up  [see  Table  III]).  However,  while  the 
radiating  mass  moved  from  370  /ig/cm  to  265  /ig/cm,  the  mass  that  was  initially  in  the  inner  plasma 
moved  from  623  /ig/cm  to  250  /ig/cm.  One  might  wonder,  therefore,  whether  there  are  advantages 
to  designing  loads  in  which  only  the  amount  of  inner,  nested,  plasma  that  ends  up  radiating  in  the 
K-shell  is  shock  heated.  Table  HI  also  shows  that  L-shell  losses  are  dependent  primarily  on  the 
amount  of  mass  in  the  outer  plasma.  Thus,  there  may  also  be  an  advantage  to  using  only  the  amount 
of  mass  in  the  outer  plasma  that  is  needed  to  shock  heat  the  radiating  core,  although  load  designs 
that  promote  emissions  from  the  outer  plasma  may  be  desired.  In  the  above  calculations,  these 
were  L-shell  titanium  emissions,  but  in  other  calculations,  perhaps  K-shell  aluminum  emissions  or 
L-shell  molybdenum  emissions  could  be  generated  in  the  outer  plasma  for  improved  simulation 
fidelity.  Finally,  further  investigation  is  needed  to  determine  the  kinds  of  load  designs  that  are 
needed  to  increase  K-shell  emissions  by  increasing  the  amount  of  core  mass  that  can  be  made  to 
participate  in  K-shell  emission.  These  investigations  have  never  been  done  for  elements  in  the 
fourth  row  of  the  periodic  table. 
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VI.  A  role  for  electron  viscosity  in  plasma  shock  heating 


I.  Introduction 

Shock  heating  of  collisional  plasmas  is  primarily  due  to  viscosity  that  dissipates 
kinetic  energy  of  the  incident  flow  directly  into  heat  in  a  so-called  viscous  sub-shock.1, 2 
The  length  scale  characteristic  of  the  viscous  sub-shock  is  of  the  order  of  the  downstream 
ion  mean  free  path,  lu ,  which  is  typically  much  less  than  (the  relevant)  hydrodynamic 

scales  of  the  surrounding  flow.  For  this  reason,  the  actual  plasma  viscosity  is  seldom,  if 
ever,  used  in  hydro  simulations  of  high-density  dynamic  plasmas.  Rather,  various  forms 
of  an  artificial  viscosity  are  commonly  used  in  Lagrangian  codes  in  lieu  of  real  viscosity. 
Conventional  justification  for  this  procedure  is  that  viscosity  is  supposed  to  affect  only 
the  shock  structure  and  its  thickness,  whereas  the  Hugoniot  conditions  that  determine  the 
shock  compression  and  the  shocked  plasma’s  energy  balance  for  a  given  shock  intensity 
are  known  to  be  independent  of  the  transport  coefficients.1'3  In  hydro  modeling,  we  are 
generally  interested  in  large-scale  phenomena;  calculations  are  not  intended  to  resolve 
structure  on  ion  mean-free-path  length  scales  (nor  they  are  generally  capable  of  resolving 
such  structures).  An  artificial  increase  in  the  shock  thickness  seems  be  a  fair  price  to  pay 
for  reducing  the  computational  effort  and  ensuring  numerical  stability  by  decreasing  the 
number  of  computational  grid  points.  This  particularly  holds  true  if  the  plasma  states 
separated  by  the  narrow  shock  front  are  accurately  modeled  when  the  true  plasma 
viscosity  is  replaced  in  the  simulation  with  some  positive-definite  artificial  expression. 

The  consistency  of  this  general  approach,  however,  comes  into  question  in  some 
situations  when  shock  waves  in  plasmas  are  modeled.  A  shock  structure  in  a  non- 
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magnetized  plasma  includes  a  long  “tongue”  of  preheat  ahead  of  the  main  viscous  jump 
in  pressure  and  velocity,  which  is  due  to  electron  thermal  conductivity.  The  length  scale 
of  this  “tongue”  is  of  the  order  of  (milme)v2lu  » /„• ,  hence  it  is  hydrodynamic  rather 

than  kinetic  scale  length.2  Since  hydro  calculations  must  account  for  and  be  able  to  model 
some  of  the  processes  occurring  within  the  shock  width,  less  arbitrariness  in  the  choice  of 
kinetic  coefficients  is  allowed.  In  particular,  because  the  electron  and  ion  temperatures  of 
the  plasma  can  be  very  different  in  shock-heated  plasmas,  separate  energy  balance 
equations  for  electrons  and  ions  must  be  used.  Then  we  have  to  decide  into  which  of 
these  two  components  of  the  plasma  should  the  kinetic  energy  dissipated  by  the  viscosity 
(artificial  or  otherwise)  be  deposited.  Does  it  all  go  to  the  ions?  Or  should  some  part  of  it 

contribute  to  the  direct  heating  of  electrons? 

For  decades,  it  was  conventionally  assumed  that  the  viscous  dissipation  heats  only 
the  ions,  the  electron  viscosity  being  of  little  consequence.  Even  the  classical  paper  by 
Braginskii4,  where  ion  viscosity  coefficients  in  fully  ionized  plasmas  are  given  for  all 
values  of  the  ion  charge  Z,  does  not  present  the  electron  viscosity  coefficients  for  Z  >  1 . 
Indeed,  for  Z  =  1  the  ratio  of  electron  to  ion  viscosity  coefficients  for  an  isothermal 
plasma,  Te=Tt,  is  very  small,  of  order  of  (me  /  m,.)1'2.  The  electron  viscosity  coefficient 

for  arbitrary  Z  calculated  in  the  same  approximation  as  used  by  Braginskii4,  was 

published  only  recently,5  correcting  the  estimate  of  the  above  viscosity  ratio  to 

(me  /  m,)1/2Z3.  Therefore,  for  Z  above  5,  the  electron  viscosity  exceeds  the  ion  viscosity 

for  equal  electron  and  ion  temperatures.  In  a  number  of  plasma  applications,  including 
intense  x-ray  generation,  shock  compression  and  heating  produce  highly  charged  ions. 
For  instance,  in  recent  experiments  with  argon6  (Z  =  18)  and  titanium7  (Z  =  22)  loads  on 
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the  “Z”  generator  at  Sandia  National  Laboratories,  substantial  quantities  of  K-shell 
radiating  helium-  and  hydrogen-like  argon  and  titanium  ions  have  been  produced. 
Advancing  further,  to  produce  K-shell  radiation  from  nickel,  copper  and  krypton  (Z 
varying  from  28  to  36),  we  have  to  heat  and  shock-compress  higher-still-Z  plasmas.  In 
these  plasmas,  the  electron  viscosity  could  considerably  exceed  the  ion  viscosity. 
Consequently,  it  must  be  taken  into  account  and  its  contribution  to  electron  heating 
estimated.  This  is  the  purpose  of  the  present  article. 

This  paper  is  structured  as  follows.  In  Section  II,  we  an  analytical  solution  of  an 
idealized  classical  shock  structure  problem  for  a  fully  ionized  high-Z  plasma,  where 
electron  viscosity  can  significantly  exceed  the  ion  viscosity.  Section  III  contains  the 
results  of  one-dimensional  (ID)  hydrodynamic  simulations  that  involve  shock 
compressions  of  high-Z  plasmas.  Here,  the  contribution  of  electron  viscosity  is 
investigated  numerically  in  a  more  realistic  time-dependent  setting.  In  Section  IV  we 
conclude  with  a  discussion  of  some  of  the  problems  in  need  of  future  investigation. 


II.  Electron  viscosity  in  the  plasma  shock  structure 

A.  Comparison  of  electron  and  ion  viscosity  coefficients 

The  classical  expression  for  the  isotropic  ion  viscosity  coefficient4  is 


Vo  ~  Voo^u  ~  Voo ' 


3m(1/27]5/2 


4 nu2  In  A e'Z 


.4-74 


(1) 


where  ,  7],  tu,  Z  and  mi  are  the  ion  number  density,  temperature  (expressed  in 
energy  units),  ion-ion  collision  time,  charge,  and  mass,  respectively.  In  A  is  the  Coulomb 
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logarithm,  and  r}^  is  a  dimensionless  coefficient  whose  value  is  given  in  Ref.  4: 
^00=0.96. 

The  electron  viscosity  coefficient  calculated  by  Whitney5  is  expressed  via  the 
similar  parameters  referring  to  the  electrons  (with  ion-ion  collision  time  tu  replaced  by 

the  electron-ion  collision  time  rfi): 


1/2t-5/2 


.  3m‘  ie 

Vo -Voo(Z)neTerei-Vco(Z)  ^nf2  In  A e4Z 


(2) 


Here,  the  dimensionless  coefficient  r\em (Z)  is  a  function  of  the  ion  charge,  Z,  given  in 


Ref.  5: 


V‘oo(Z)  = 


1.81Z(Z2  +2.82Z  +  1.343) 
Z3  +  4.434Z2  +5.534Z  +  1.78 


(3) 


For  Z  =  1,  we  obtain:  r\eM  =  0.733 ,  the  value  calculated  by  Braginskii.4  From  (l)-(3),  we 


find  the  ratio  of  viscosity  coefficients: 


ai= 

Vo 


M 

1/2 

(t  ^ 

*  e 

mu 

\  u  ) 

'  J 

^(Z)Z3 


(4) 


<(2A)1/2  ’ 

where  mu  is  the  atomic  mass  unit  (mu  =  mproton ),  A  is  the  ion  atomic  weight.  For  a  fully 
ionized  plasma  with  not-too-high  atomic  number,  one  can  approximate  A  =  2Z ,  to 
obtain1^ 


Vi  2 


V'VZ7^ 

T, 


\mu  J 


5/2 


fse\5'2 


.  1.17X10-'  * 


\  1  J 


vSo(Z) 

V‘oo 


(5) 


f  in  Ref.  5,  where  this  ratio  was  first  presented,  a  factor  of  4/3  was  inadvertently  omitted 
from  the  expression  for  the  ion  pressure  tensor.  Here  we  present  the  corrected  formula  for 
the  electron-to-ion  viscosity  ratio. 
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where  and  if  are  thermal  energy  densities  of  the  plasma  components.  To  see  why  the 

ratio  of  viscosity  coefficient  depends  on  Z  in  this  way,  one  can  consult  elementary  kinetic 
theory.4  Recall  that  the  viscosity  coefficient  for  any  sort  of  particle  is  estimated  to  be  of 
order  nTx  =  nT In  vTo ,  where  n  is  the  colliding  particle  number  density, 

vT  =  (3 T /m)u2  is  the  average  thermal  velocity,  n  is  the  number  density  of  the  particles, 
against  which  the  collisions  occur.  The  collision  time  r  is  defined  in  terms  of  o ,  the 
average  cross  section  of  such  collisions.  For  charged  particles,  c  is  proportional  to  the 
product  of  both  charges  squared  times  the  Coulomb  logarithm  divided  by  the  temperature 
squared.  In  high-Z  plasmas,  there  are  Z  more  electrons  than  ions,  and  the  electron  mean 
free  path  is  longer  than  the  ion  mean  free  path.  Indeed,  for  ion-ion  collisions,  n  =  n  =nn 

and  <JU  °c  Z4 ,  which  leads  to  (1).  The  electron-ion  collision  cross  section  aei  is  less  than 

<t(7  by  a  factor  of  Z2 ;  additionally,  electrons  outnumber  ions  by  a  factor  of  Z 

(n  =  ne  =  Znn  n  =  ni ).  The  factor  Z3  in  (4)  comes  from  this  dependence  on  n,  n,  and  o. 

Combining  this  factor  with  the  factor  A~U1  oc  Z~l/2  coming  from  the  ratio  of  thermal 
velocities,  we  arrive  at  Eq.  (5).  In  high-Z  plasmas,  most  of  the  thermal  energy  typically  is 
in  the  electron  component,  and  we  obtain  from  (5)  that  electron  viscosity  exceeds  ion 
viscosity  when  %l%t>  5 .  The  ratio  of  viscosity  coefficients  vs.  Z  given  by  (4)-(5)  is 

plotted  on  Fig.  1  for  fully  ionized  plasmas  with  Te  =Tr  We  see  that  r)eQ  >rj‘0  for  Z  >  5 

(carbon  and  higher).  The  ratio  rjg  lr\‘0  exceeds  Z,  implying  greater  viscous  heating  per 

electron  than  per  ion,  for  Z  >  13  (silicon  and  higher).  For  argon  and  higher-Z  plasmas 
stripped  to  helium-  and  hydrogen-like  states,  as  in  experiments  of  Refs.  6,7,  the  electron 
viscosity  is  substantially  greater  than  the  ion  viscosity. 
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B.  Classical  plasma  shock  structure  with  electron  viscosity  included 

Let  us  find  out  how  the  presence  of  substantial  electron  viscosity  affects  plasma 
shock  heating.  For  this  purpose,  we  modify  the  well-known  solution  of  a  classical 
collisional  plasma  shock  structure  problem,  originally  studied  in  Refs.  8-10  (see  also 
Section  2.1  of  Ref.  2),  to  include  the  electron  viscosity.  In  this  Section,  we  consider  a 
fully  ionized  plasma  with  constant  ion  charge  Z,  ignoring  the  ionization  kinetics  and 
thereby  excluding  contributions  of  ionization  and  excitation  from  the  energy  balance,  as 
well  as  radiation  effects,  which  are  all  significant  in  high-Z  non-fully-ionized  plasmas. 
These  are  accounted  for  in  the  simulations  that  are  reported  in  Section  III.  Additionally, 
we  neglect  the  ion  thermal  conduction  whose  characteristic  length  scale  in  the  shock 
structure  we  do  not  attempt  to  resolve.  Here  we  limit  ourselves  to  an  idealized  problem  to 
single  out  and  highlight  specifically  the  magnitude  of  the  effect  of  electron  viscosity. 

To  find  a  steady-state  shock  structure,  we  conventionally  pass  to  the  reference 
frame  where  the  shock  front  is  stationary  and  write  down  the  steady-state  ID  fluid 
equations,4  which  in  this  case  simply  reduce  to  integrated  form  of  conservation  equations 

l  2 

of  mass,  momentum,  and  energy  ' 

(6) 

nu  =  «,«, ,  v  1 

/wijH2  +n(ZTe  +Tj)——(rjo  +Tl'o)—^  ~  n\miul  +  (1+Z)n,7j ,  (2) 

+-nu(ZTe  +  Ti)-Ke^-%e0 

2  2  dx  3  dx  (8) 

=  ■in1m(M13  +  -|(1  +  Z)n{uJ{ 

L  +* 

and  the  heat  balance  equations  for  the  ions 
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(9) 


3  dT.  du  4  /  du) 

-nu — L  +  nT  —  =  -V0  — 

2  dx  dx  3  ydx  j 


+Zn3J^LzL 

m.  r„ 


Here  n  is  the  ion  number  density  (due  to  plasma  quasi-neutrality,  the  electron  number 
density  ne  -Zn ,  and  the  velocity  u  is  the  same  for  both  electrons  and  ions); 


Ke  =  Kl(Z)neTex e  / me  is  the  electron  heat  conduction  coefficient 


(re  =  3my2r,3/2 /4a/27T  In A<?4Z2n ,  and  the  dimensionless  coefficients  K*(Z)  are 

tabulated  in  Ref.  4.)  Subscript  1  refers  to  the  values  of  the  flow  variables  in  the  pre-shock 
uniform  flow  far  upstream  from  the  shock  front,  where  all  the  gradients  vanish;  in 
particular,  w,  is  the  shock  speed  defined  with  respect  to  the  plasma  ahead  of  it.  But  the 
flow  variables  also  vanish  in  the  post-shock  state,  far  downstream,  which  we  would 
designate  by  subscript  2.  The  right-hand  sides  of  Eqs.  (6)-(8)  must  have  the  same  values 
whether  they  are  expressed  via  the  parameters  characteristic  of  the  pre-shock  (subscript 
1)  or  the  post-shock  (subscript  2)  state.  From  this  observation,  we  derive  the  relations 
between  the  pre-  and  post-shock  state  parameters: 


u2  _  nx  _  A/,2  +3 
~  n2~  4  M2  ~  2’ 

T,  (M2+3)(5M,2-1)  _ 
T]~  16  M,2 


(ID 


(12) 


where  Mx  =  m,  •  [5(1  +  Z)TX  I3mi  ]~l/2  is  the  Mach  number  characteristic  of  the  shock  wave. 
These  are,  of  course,  conventional  Hugoniot  conditions  for  an  ideal  gas  with  y  =  5/3 . 
Since  Mx  >  1 ,  we  have  U2<  1  and  ©2  >  1 :  the  shock  wave  compresses  and  heats  the 
plasma. 

To  obtain  the  shock  structure  solution,  it  is  convenient  to  use  normalized  variables 
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(13) 


Then  Eq.  (6)  is  satisfied  automatically  if  nx  In  equals  U  and  Eqs.  (7)-(9)  are  reduced  to 


U-1+- 


5M, 


ze,+e,.  t 

(1  +  Z)U 


a  Jf/  n 
-A  —  =  0, 

dx 


-(U2-  1)  +  -^T 

2  2  M,2 


Z0  +0,  ^ 


Z  +  l 


-1 


A  T7dU  .  d@ 

-AvU- - Ar  — ^  =  0, 

dx  dx 


(14) 

(15) 


3rf0L+0Lrff/  =  5O  +  ZKM2A 

2&  f/  dx:  3  r?o+t7o 


f— 1 

,  ^  > 


(16) 


where  the  length  scale  functions  Av,ATe,  and  Ar  of  the  variables  U,  Qe,  0,  (and 
parameters  Mx,  Z),  which  are  characteristic  of  viscosity,  electron  thermal  conductivity 
and  temperature  relaxation,  respectively,  are  defined  as 

A  A  JOL  A  t  £/.  (17) 

*  ’  Te  3  r  o  ^7 _  1  e 


3Zm . 


3  m^w. 

We  are  interested  in  a  solution  of  Eqs.  (14)-(16)  with  the  Hugoniot  boundary  conditions: 
U,  ©,,  ©,  ->1  at  U->U2,  0e,  0,  ->©2  at  x ->+<*>.  (18) 

The  upstream  and  downstream  points  at  infinity,  where  all  the  derivatives  vanish, 
correspond  to  singular  points  of  the  autonomous  dynamic  system  (14)-(16).  Therefore  the 
solution  sought  for  represents  a  singular  trajectory  that  joins  the  two  singular  points  (18). 

The  solution  is  obtained  using  the  multiple-scale  method,  which  is  described  in 
details  in  Refs.  2  and  10.  The  largest  scale  length  turns  out  to  be  ATe ,  which  is  one 


associated  with  electron  thermal  conduction,  and  which  is  comparable  to  the  temperature 
relaxation  scale,  Af .  Consequently,  by  assuming  that  all  the  derivatives  dldx,  are 
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~  1/ ATe ,  by  neglecting  the  small  terms  of  order  of  Av  /ATe ,  and  introducing  a 


dimensionless  length  variable, 

>  £=*/Ar<,,, 

one  can  reduce  the  system  (14)-(16)  to  an  algebraic  relation  for  0, 


=  -Z0e+(Z  +  l) 

©  +—M?U2  - 

|m2+iV 

3 

U  )  J 

and  to  the  following  two  equations  for  U  and  ©e : 


(19) 


(20) 


dU 


© 

c 

( 5  2  ) 

3Z(U -l)(U2 -U)  +  s—j 

0  +-m2u2- 

— M2  +1 

U 

^  T  J  * 

3 

l3  J 

-Z%  +  -(Z  +  l)(-8M2f/+5M2+3) 
U2  6  1  1 


(21) 


^  =  2©;5/2(t/-l)(C/,-f7),  (22) 

where  s  =  27kt<0(Z)Z2 /25(1 +  Z)2M,4 . 

An  asymptotic  analysis  of  the  system  (21)-(22)  in  the  vicinity  of  its  singular 
points,  1  and  2'  [see  (18)]  reveals  that  the  singular  point  1  is  a  saddle  points  for  any 
M,  >  1 .  Consequently,  a  unique  phase  trajectory  on  the  phase  plane  ( U ,  ©e)  approaches 
this  singular  point  (1,  1)  at  £  — »  from  the  domain  of  higher  temperatures  and 

densities.  A  shock  structure  solution  would  exists  and  be  physically  meaningful  (stable, 
evolutionary2)  if  the  same  phase  trajectory  can  be  found  to  approach  to  the  singular  point 
2  (U2 , 02)  at  £  — » +°° .  This  is  the  case  when  the  singular  point  2  is  a  node  attracting  at 
£  —>  +°o  all  the  trajectories  that  pass  in  its  vicinity.  These  latter  conditions  holds  for 
weak  shock  waves,  defined  by  the  inequality 
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M,  <  Mcr 


9Z  +  10 
5(Z  +  2) 


(23) 


The  structure  of  a  weak  shock  wave  is  determined  by  electron  thermal  conduction, 
viscosity  playing  no  role  at  all.  Expression  (23)  for  a  critical  Mach  number  is  a  simple 
generalization  of  Shafranov’s  result9  (which  was  derived  for  Z—  1,  implying 
Mcr  =  -s/19/15  =  1.125 )  to  arbitrary  Z.  Weak  shock  profiles  for  high-Z  plasmas  are  no 

different  from  those  discussed  in  Refs.  9-1 1,  2  (in  particular,  electron  and  ion 
temperatures  are  very  close  throughout  the  shock),  and  we  will  not  describe  them  any 
further  here.  Rather,  we  are  interested  in  the  strong-shock  case,  where  electron  viscosity 
does  make  a  difference. 

For  Mx  >  Mcr(Z) ,  the  downstream  singular  point  2  is  a  saddle  point,  too.  In  other 
words,  only  one  phase  trajectory  approaches  it  as  £  ->  +°°  from  the  domain  of  lower 
temperatures  and  densities,  and  this  is  not  the  same  trajectory  that  approaches  the 
downstream  singular  point  1  in  the  limit  £  — >  .  As  a  consequence,  the  shock  structure 

defined  by  (20)-(22)  contains  a  “viscous  sub-shock”,  which,  to  zeroth  order  in  Av  /  Are ,  is 
simply  a  zero-thickness  discontinuity,  a  jump  from  the  former  trajectory  to  the  latter.  If 
the  length  scale  of  the  jump  is  Av,  then  d/dx  ~  1/  A„ .  In  this  case,  Eq.  (15)  ensures  that 
electron  temperature  in  the  viscous  sub-shock  is  essentially  constant9, 10  due  to  the 


smallness  of  the  ratio  A„  /  Are . 

To  obtain  the  Hugoniot  conditions  for  the  viscous  sub-shock,  one  can  neglect  the 
thermal  relaxation  term  in  the  right-hand  side  of  Eq.  (16)  (its  order  is  determined  by  the 
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small  ratio  Av  /  Ar).  Then,  assuming  constant  electron  temperature  0e  =  0M  and 
simplifying  this  equation  with  the  aid  of  (14),  we  reduce  it  to 


dU  3(05/2  + 


^L(Z+WJ^L+1\Z+1)+Z^ 

3  3  U  U 


where  n  denotes  the  right-hand  side  of  Eqs.  (4)  or  (5)  at  Te  =  Tr  Equation  (24)  is  to  be 
integrated,  from  the  point  with  the  initial  condition  0,  =  ©ifl  at  U  =  U a  to  the  point 
U  -  Ub,  0  =  Q  b  where  the  velocity  derivative  in  the  left-hand  side  of  Eq.  (14)  vanishes. 


At  this  latter  point. 


Ub-l  +  - 


3  ZSea  +  Qib 


5M,  (1  +  Z)Ub 


-1=0. 


Points  a  and  b  are  thus  related  by  the  Hugoniot  condition  for  the  viscous  sub-shock. 

In  an  analysis  in  which  electron  viscosity  is  neglected  by  setting  n  -  0 ,  one  can 
explicitly  integrate  Eq.  (24).  Then,  by  substituting  the  result  of  the  integration  into  (25), 
one  can  reduce  the  subshock  Hugoniot  condition  to  a  single  transcendental  equation:10 


(Ua-Ub)(l+U2-Ua-U„)- 


iom;(z+i)  u„ 


The  structure  of  a  shock  with  electron  viscosity  included  is  calculated  similarly  to 
Refs.  9,  10.  We  start  in  the  vicinity  of  singular  point  1  and,  on  the  phase  plane  (U,  0e) , 
plot  the  only  trajectory  of  the  system  (21)-(22)  that  is  outgoing  from  point  1  at  £  — >  to 

the  domain  of  higher  temperatures  and  pressures.  This  part  of  the  solution  represents  the 
electron  pre-heat  zone  ahead  of  the  viscous  shock.  For  each  point  a  on  this  trajectory,  we 


find  the  corresponding  point  b,  which  is  connected  to  it  by  the  Hugoniot  relation  for  the 
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viscous  sub-shock,  as  explained  above.  As  the  point  a  moves  along  the  pre-heat 
trajectory,  it  maps  out  a  locus  of  the  corresponding  points  b  on  the  phase  plane  (Fig.  2). 
The  shock  structure  equations  permit  a  viscous  subshock  transition  (occurring  at  constant 
electron  temperature,  i.e.  a  horizontal  segment  on  Fig.  2)  from  any  point  a  to  each 
corresponding  point  b.  However,  only  one  b  point  lies  on  the  solution  trajectory.  Starting 
from  the  correct  point  b,  the  solution  to  Eqs.  (20)-(22)  must  approach  the  downstream 
singular  point  2  at  £  ->  +°° ,  and  this  can  happen  only  along  a  single  unique  phase 
trajectory.  To  find  the  only  trajectory  that  enters  point  2  at  £  +°° ,  we  start  in  the 

vicinity  of  point  2,  and  integrate  Eqs.  (21)-(22)  backward  over  £ .  After  locating  the 
intersection  of  this  trajectory  on  the  phase  plane  (U,  0e)  with  the  locus  of  points  b,  we 
can  construct  the  whole  solution  1  a  ->  b  -»  2 ,  which  satisfies  the  general  Hugomot 
conditions  (13). 

The  Hugoniot  conditions  for  the  viscous  sub-shock  are  different,  depending  on 
whether  or  not  the  electron  viscosity  is  accounted  for.  Therefore,  for  the  given  outgoing 
trajectory  1 4  a ,  we  can  plot  plane  two  curves  on  the  phase,  each  representing  a  locus  of 
points  b :  one  for  the  case  when  the  electron  viscosity  is  accounted  for,  and  the  other  for 
the  case  when  it  is  ignored  (dotted  and  dashed  locus  curves  on  Fig.  2,  respectively.)  From 
each  point  ax  there  is  a  possible  viscous  shock  transition  (represented  by  the  dotted 
horizontal  segment)  either  to  bx  or  to  b2 ,  depending  on  whether  or  not  the  electron 
viscosity  is  accounted  for.  Each  of  the  two  locus  b  curves  intersects  with  the  only 
trajectory  entering  2  at  £  -» +~  at  a  different  location.  Note  that  solutions  of  Eqs.  (20)- 
(22)  in  the  vicinities  of  pre-  and  post-shock  states  are  exactly  the  same,  no  matter  whether 
or  not  the  electron  viscosity  is  accounted  for.  The  only  difference  between  them  is  in  the 
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location  and  strength  of  the  viscous  sub-shock.  Not  surprisingly,  with  electron  viscosity 
included,  the  pre-heat  zone  is  a  little  longer,  and  the  electrons  are  pre-heated  to  a  higher 
temperature,  whereas  the  temperature  relaxation  zone  is  accordingly  shorter,  and  the  ions 
in  the  viscous  sub-shock  are  less  over-heated  compared  to  the  electrons.  For  the  example 
of  Fig.  2  (Z  =  22,  M,  =  5>Mcr(22)  =  1.32)  we  find:  without  electron  viscosity  the  peak 
ion  temperature  in  the  viscous  sub-shock  exceeds  the  peak  (that  is,  post-shock)  electron 
temperature  by  the  factor  3.8,  whereas  once  electron  viscosity  is  taken  into  account  this 
ratio  is  reduced  down  to  1.7.  The  shock  profiles,  for  both  cases,  are  plotted  in  Fig.  3. 

Why  are  the  ions  still  overheated  when  the  viscous  energy  release  per  electron  is 
higher  than  per  ion  (by  a  factor  of  ~5,  in  this  example)?  The  reason  is  that  thermal  energy 
is  drained  from  the  electrons  via  thermal  conduction  upstream.  The  heat  released  in  the 
viscous  sub-shock  is  instantly  spent  to  heat  the  electrons  upstream.  On  the  contrary,  the 
length  scale  characteristic  of  the  ion  heat  conduction  (which  is  ignored  here)  is  of  the 
same  order  as  the  thickness  of  the  viscous  sub-shock.  Thus,  although  the  viscosity 
generates  less  heat  per  particle  in  the  ion  component  of  the  plasma,  the  ions  keep  this 
energy,  only  slowly  sharing  it  with  electrons  behind  the  viscous  sub-shock,  in  the 
temperature  relaxation  zone.  Figure  4  shows  how  the  ion  overheating  (characterized,  as 
above,  by  the  ratio  of  the  peak  ion  temperature  to  the  post-shock  plasma  temperature)  for 
a  strong  shock  wave  (M,  =  10 )  varies  as  Z  is  varied  from  13  (aluminum)  to  36  (krypton). 
If  the  electron  viscosity  is  ignored,  then  overheating  is  seen  to  be  a  linearly  increasing 
function  of  Z  for  a  given  Mach  number.  With  electron  viscosity  taken  into  account, 
overheating  is  virtually  constant,  about  1.7,  in  the  range  of  Z  shown.  The  physical 
meaning  of  both  scalings  is  quite  clear.  If  electron  viscosity  is  neglected,  then  each  ion 
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dissipates,  in  the  viscous  sub-shock,  its  own  kinetic  energy,  which  is  proportional  to  the 
ion  mass  for  a  given  shock  velocity.  Hence,  the  dissipated  energy  per  ion  scales  as  A,  or, 
equally,  as  Z;  whereas,  far  downstream,  this  energy  is  shared  with  Z  electrons,  thus 
reducing  the  ion  temperature  by  a  factor  of  Z.  If  electron  viscosity  is  included  in  the 
energy  dissipation,  this  energy  sharing  between  an  ion  and  its  Z  electrons  starts  already  in 
the  viscous  sub-shock,  making  the  overheating  virtually  independent  of  Z. 

III.  Effect  of  electron  viscosity  on  radiating  Z-pinch  plasmas 

Let  us  study  in  a  more  realistic  and  time-dependent  hydrodynamic  calculation  one 
example  of  a  physical  situation  involving  plasma  shock  waves  where  taking  the  account 
of  the  electron  viscosity  makes  a  difference.  Shock  heating  is  known  to  be  a  major 
mechanism  for  converting  kinetic  energy  of  imploding  Z-pinch  plasmas  into  thermal 
energy,  making  the  kinetic  energy  available  for  producing  hard  x-ray  radiation.  To 
obtain  K-shell  emission  in  the  multi-keV  range,  the  ions  have  to  be  stripped  down  to 
high-Z  helium-  or  hydrogen-like  states,  e.  g.,  we  need  to  reach  Z  =  20  -  21  or  34  -  35  in 
order  to  produce  either  titanium  or  krypton  K-shell  radiation,  respectively.  The  Z-pinch 
plasma  is  shocked  at  least  once,  in  the  case  of  an  empty  annular  shell  imploding  onto  the 
axis.  In  the  case  of  a  gas-puff  load,  whether  it  be  a  gaseous  “shell”  or  a  uniform  fill,  there 
is  always  some  gas  near  the  axis,  so  a  converging  shock  wave  reaches  the  axis  ahead  of 
the  current  sheath,  after  which  the  plasma  is  re-shocked  in  a  reflected  wave.  In  wire-array 
implosions,  the  so-called  precursor  plasma  also  reaches  the  axis  prior  to  the  main  plasma 
mass,14, 15  hence  plasma  near  the  axis  can  be  shocked  twice,  too.  Three  or  more  plasma 
shock  fronts  could  be  produced  with  more  advanced  loads,  like  double  gas  puffs16  or 
nested  wire  arrays.17 
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With  ionization  starting  in  the  pre-heat  zone,2  even  as  the  first  (and  particularly 
the  subsequent)  viscous  sub-shocks  compresses  and  heats  a  high-Z  plasma,  electron 
viscosity  becomes  important.  Its  role  is  illustrated  by  the  results  of  one-dimensional  (1-D) 
simulations  of  titanium  and  krypton  Z-pinch  implosions  given  below.  We  did  not 
undertake  a  systematic  and  thorough  study  of  how  electron  viscosity  affects  Z-pinch 
implosions.  Rather,  our  present  results  represent  only  a  first  look  at  how  the  electron 
viscosity  influences  K-shell  emission,  and  about  how  this  influence  in  certain  implosions 
is  affected  by  the  atomic  number  of  the  plasma. 

We  simulated  titanium  and  krypton  load  implosions  with  the  1-D  radiative-MHD 
code  described  in  Ref.  18.  The  calculations  were  carried  out  using  from  60  to  100 
Lagrangian  zones.  A  collisional-radiative  equilibrium  (CRE)  ionization  dynamics  was 
employed  in  all  cases.  The  standard  form  of  artificial  viscosity  was  used.  The  heating  it 
produced  was  either  conventionally  put  entirely  into  the  ions  or  divided  between 
electrons  and  ions  according  to  the  ratio  of  their  physical  viscosity  coefficients,  Eq.  (4), 
and  the  results  of  the  two  procedures  were  compared. 

Our  calculations  were  done  for  two  kinds  of  imploding  loads:  1)  uniform  fills  of 
krypton  and  titanium,  and  2)  nested  titanium  shell-on-fill  plasmas.  In  the  first  set  of 
calculations,  the  initially  uniform  density  fills  of  krypton  and  titanium  plasmas  were  at  a 
temperature  of  4  eV.  The  mass  of  the  fill,  its  initial  density,  and  the  energy  input  were 
varied  so  that  comparisons  of  the  ion  temperature  profiles  and  K-shell  yields  and  powers 
could  be  made  with  and  without  electron  viscosity  in  the  calculations.  A  nested  titanium 
plasma  geometry  was  set  up  by  moving  an  outer  portion  of  the  above  described  fill 
farther  out  from  the  inner  portion,  and  creating  a  small  mass  of  reduced  density  (by  a 
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factor  of  10)  between  the  inner  and  the  outer  plasmas.  The  loads  were  driven  by  a  linearly 
rising  current  pulse  that  was  terminated  when  a  prescribed  energy  input  was  achieved. 

Fill  masses  were  varied  between  1.25  and  10  mg/cm.  Most  of  this  interval  physically 
corresponds  to  next  generation  of  multi-MA  current  drivers  capable  of  producing  krypton 
K-shell  emission  (energy  delivered  to  the  plasma  up  to  ~12  MJ/cm,  current  up  to  ~80 
MA).  The  nested  calculations  were  performed  with  a  plasma  mass  of  1.25  mg/cm,  which 
corresponds  to  an  amount  of  plasma  mass  that  can  be  successfully  imploded  on  the  Z 
accelerator  to  velocities  that  are  needed  to  effectively  excite  the  K-shell  of  titanium 

7 

(energy  delivered  to  the  plasma  ~0.5  MJ/cm,  current  ~18  MA.) 

When  electron  viscosity  is  ignored  and  all  of  the  viscous  heating  is  put  into  ions, 
the  rate  of  ion  heating  exceeds  the  electron  rate,  which,  in  this  case,  is  due  to  adiabatic 
compression  heating  alone.  Thus,  the  ions  become  significantly  hotter  than  the  electrons 
(cf.  Fig.  3,  dashed  lines),  in  all  the  classical  shock  structure  calculations  that  include  ion 
viscosity  only.  As  the  plasma  become  highly  ionized,  the  addition  of  electron  viscosity 
should,  therefore,  reduce  ion  temperatures  (Fig.  3,  solid  lines)  and  raise  electron 
temperatures  slightly,  which  suggests  the  possibility  of  converting  some  of  this  energy 
into  extra  K-shell  emission.  However,  this  possibility  depends  on  the  plasma  conditions 
that  are  generated  during  the  actual  dynamics  of  a  time-evolving  implosion,  since  the 
viscosity  coefficient  ratio  is  substantially  reduced  when  7]  is  allowed  to  become  much 

larger  than  Te  in  the  early  phases  of  the  implosion  Moreover,  in  a  dynamic  implosion, 

radiation  cooling  becomes  a  significant  drag  on  the  electron  temperature  that  acts  on  the 
ions  only  indirectly  through  electron-ion  collisions.  Thus,  the  ions  can  get  hotter  than  the 
electrons  when  radiation  losses  are  high  and  the  electron-ion  coupling  is  too  weak.  Under 
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conditions  of  too  low  an  ion  density,  and  too  small  an  electron-ion  collision  rate,  the 
rerouting  of  ion  thermal  energy  into  K-shell  emission  is  inhibited. 

Figure  5,  6  show  our  results  for  a  10  mg/cm  Kr  uniform  fill  with  an  initial  4  cm 

radius  driven  by  an  energy  input  ~12  MJ/cm.  This  corresponds  to  7]  =2,  where  tj 
stands  for  the  energy  delivered  per  ion  divided  by  the  energy  required  to  strip  the  ion  into 
He-like  excited  states.  Figure  5  demonstrates  the  electron  and  ion  temperature  profiles 
during  the  run-in  phase,  at  80  ns  (a)  and  near  the  peak  of  the  K-shell  emission,  at  113.5  ns 
(b).  We  see  that  the  inclusion  of  electron  viscosity  produces  significant  reductions  in  the 
ion  temperature.  Reduced  ion  heating  and  the  corresponding  increase  in  electron  heating 
led  to  a  20%  increase  in  the  K-shell  energy  output,  making  the  total  K-yield  approach  4 
MJ/cm,  see  Fig.  6.  All  of  the  K-shell  emission  originated  from  the  core  of  the  compressed 
plasma.  The  outer  portion  of  plasma,  which  contains  a  sizable  fraction  of  the  plasma 
mass,  emitted  only  in  the  L-shell.  Prior  to  the  K-shell  emission  [Fig.  5(a)],  the  electron 
temperature  profile  shows  only  slight  changes  due  to  electron  viscosity  in  contrast  to  the 
ion  temperature  profile  since  there  are  Z  times  as  many  electrons  to  share  the  energy  as 
there  are  ions. 

For  a  10  mg/cm  uniform-fill  of  titanium  imploded  with  an  initial  radius  of  2  cm, 
similar,  but  smaller  reductions  in  ion  temperature  were  seen,  as  illustrated  by  Figs.  7,  8. 
Figure  7  shows  the  ion  temperature  profiles  near  the  peak  of  compression  for  77*  =1.8, 
which  corresponds  to  the  current  ~40  MA  and  the  energy  delivered  to  the  pinch  about  3 
MJ/cm,  almost  30%  of  which  is  emitted  in  the  K-shell.  The  reductions  in  the  ion  thermal 
energy  again  showed  up  as  increases  in  K-shell  emission,  although  to  a  somewhat  smaller 
extent,  8%  in  Fig.  8.  Doubling  the  energy  input  twice,  to  rf  =  3.6 ,  and  then  to  rf  =  7.2 
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(this  was  achieved  by  reducing  the  rise  time  of  the  current),  we  find  corresponding 
increases  in  the  K-shell  emissions  of  the  same  order,  15%  and  10%,  respectively. 

Figure  9  shows  the  ion  temperature  profiles  obtained  for  implosion  of  a  nested 
titanium  wire  array.  The  nesting  used  in  the  calculations  was  done  by  dividing  the  1.25 
mg/cm  uniform-fill  plasma  into  two  parts  and  by  moving  the  outer  part  in  radius,  by 
roughly  a  factor  of  two.  Since  one  cannot  have  a  perfect  vacuum  between  plasmas  in  a 
simulation,  the  calculations  were  zoned  so  that  the  density  between  the  inner  and  the 
outer  plasmas  could  be  reduced  by  a  factor  of  10.  Therefore,  the  calculations  reported 
here  are  of  the  shell-on-fill  rather  than  shell-on-shell  variety,  thereby  representing  a  very 
rough  model  of  the  actual  “Z”  experiments.  Still,  even  though  a  perfect  vacuum  did  not 
exist  between  the  two  plasmas  in  the  calculation,  the  collision  time  is  well  defined.  In  the 
run  illustrated  by  Fig.  9,  the  outer-to-inner  mass  ratio  was  taken  to  be  2,  the  outer  plasma 
radius  4.2  cm,  energy  input  0.42  MJ/cm  (which  translates  into  rf  =  2  and  current  of  16 
MA),  total  number  of  zones  100.  The  temperature  profiles  are  shown  in  Fig.  9  at  the 
onset  of  the  K-shell  emission.  Here  again,  the  inclusion  of  electron  viscosity,  by  reducing 
the  ion  temperature,  allowed  an  additional  20  kJ  of  K-shell  x  rays  to  be  extracted  from 
the  plasma,  which  increased  the  K-shell  yield  by  about  10%. 

IV.  Conclusions 

Two  methods  were  used  to  evaluate  the  influence  of  electron  viscosity  on  the 
shock  heating  of  highly  ionizable  moderate  atomic  number  plasmas.  Classical,  analytic 
shock  structure  calculations  were  first  carried  out  to  quantify  the  amount  by  which  the  ion 
temperature  discontinuity  across  the  unresolved  sub-shock  structure  of  a  plasma  shock 
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wave  was  reduced  by  taking  electron  viscosity  into  account.  The  strength  of  this 
calculation  is  that  it  can  make  use  of  real  classical  viscosities  and  accurately  resolve  the 
shock  structure  that  is  determined  mainly  by  the  electron  thermal  conduction.  On  the 
other  hand,  its  weaknesses  are  that  it  cannot  treat  convergent  shock  waves  and  that  it 
must  assume  a  steady  flow  and  a  fully  ionized  plasma  on  both  sides  of  the  shock  front.  A 
second  method  involved  the  use  of  computer  calculations  to  quantify  how  electron 
viscosity  influences  ion  heating  in  strongly  driven  and  shock-heated  Z-pinch  plasmas. 

The  strengths  of  these  calculations  were  that  they  could  treat  both  time-evolving,  ionizing 
and  radiating  plasmas,  and  convergent  shock  wave  dynamics.  However,  their  weaknesses 
were  that  they  could  not  resolve  the  shock  structure,  nor  could  they  properly  treat  the 
non-equilibrium  ionization  dynamics  across  the  sharp  temperature  and  density  gradients. 
Consequently,  artificial  viscosities  were  used  in  these  calculations,  which  do  not  include 
any  of  the  Z  or  temperature  dependences  of  the  real  viscosities.  These  dependences  were 
included  though  in  the  way  the  artificial  viscosity  heating  was  divided  between  electrons 
and  ions. 

We  have  demonstrated  that  the  inclusion  of  previously  neglected  electron 
viscosity  can  significantly  affect  the  results  of  hydro  simulations  of  shock-heated  high-Z 
plasmas.  The  main  effect  is  a  considerable  reduction  of  the  ion  over-heating  in  the  shock 
waves.  Hence  we  predict  that,  for  instance,  direct  measurements  of  ion  temperature  in 
high-current  Z-pinches  (e.  g.,  those  based  on  Doppler  broadening)  are  likely  to  show  the 
ion  temperatures  that  are  much  lower  than  those  expected  from  simulations  in  which 
electron  viscosity  is  not  taken  into  account.  A  more  subtle  effect  is  the  corresponding 
rather  small  increase  in  the  electron  temperature.  If  electron  viscosity  is  neglected,  the  ion 
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overheating  in  the  viscous  sub-shock  scales  approximately  as  Z  (Fig.  4),  the  same  as  the 
number  of  electrons  per  ion.  By  including  the  electron  viscosity,  we  split  most  of  the 
ion’s  extra  thermal  energy  between  all  Z  electrons.  Consequently,  when  electron  viscosity 
is  added,  the  resulting  relative  increase  of  the  electron  energy  (temperature)  must  be  a 
relatively  weak  function  of  Z.  This  conclusion  is  consistent  with  our  simulation  results, 
where,  in  the  range  of  Z  between  22  and  36,  increases  of  the  K-shell  output  due  to  the 
extra  energy  diverted  by  the  electron  viscosity  was  found  to  vary  only  between  10  and 
20%.  Therefore  we  can  conclude  that,  whenever  a  shock  compression  of  high-Z  plasmas 
is  modeled,  the  electron  viscosity  should  be  taken  into  account.  This  is  particularly  true  if 
we  are  interested  in  more  accurately  predicting  the  ion  temperature,  or  if  our  results  are 
sensitive  to  the  value  of  the  electron  temperature. 
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Figure  captions 

Fig.  1 .  Ratio  of  electron  and  ion  viscosity  coefficients  atTe=Tr  Points  are  plotted  for  the 


exact  values  of  atomic  masses  A(Z) ,  curves  correspond  to  the  approximation  A  =  2 Z. 

For  Z  =  1,  we  take  A  =  2  (deuterium). 

Fig.  2.  Phase  plane  (U,QJ .  Phase  trajectories  leaving  the  upstream  point  1  and  entering 

the  downstream  point  2  correspond  to  the  electron  pre-heat  zone  and  to  the  electron-ion 
temperature  relaxation  zone,  respectively.  Horizontal  segments  are  viscous  sub-shocks 
permitted  by  the  Hugoniot  conditions,  dotted  and  dashed  curves  -  loci  of  points  b  with 
and  without  electron  viscosity,  respectively. 

Fig.  3.  Shock  profiles  for  a  =5  strong  shock  wave  in  a  fully  ionized  titanium  plasma, 
Z  =  22.  Here,  normalized  density  equals  \IU ,  see  (13). 

Fig.  4.  Ion  overheating  in  a  M,  =  10  strong  shock  wave  in  a  fully  ionized  high-Z  plasma, 
as  a  function  of  Z. 

Fig.  5.  Ion  temperature  profiles  for  a  Kr  uniform-fill  implosion  calculated  without  (A) 
and  with  electron  viscosity  (B).  (a)  during  the  run-in  phase;  (b)  near  the  peak  of  K-shell 
emission. 

Fig.  6.  Time  history  of  Kr  K-shell  yield  for  the  conditions  of  Fig.  5 

Fig.  7.  Ion  temperature  profiles  for  a  Ti  uniform-fill  implosion  calculated  without  (A)  and 

with  electron  viscosity  (B).  (a)  during  the  run-in  phase;  (b)  near  the  peak  of  K-shell 

emission. 

Fig.  8.  Time  history  of  Ti  K-shell  yield  for  the  conditions  of  Fig.  7. 

Fig.  9.  Ion  temperature  profiles  for  a  Ti  uniform-fill  implosion  calculated  without  (A)  and 
with  electron  viscosity  (B)  near  the  peak  of  K-shell  emission. 
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VII.  Improved  Bremsstrahlung  Radiation  Production  by 
Electrode  Decontamination 


The  focus  of  this  work  is  on  the  need  to  demonstrate  on  Decade  the  improvements  seen  in 
bremsstrahlung  radiation  production  from  heating  the  tantalum  anode  to  white-hot 
temperatures.  Last  year  a  test  on  the  Gamble  II  machine  at  NRL,  white-hot  tantalum 
anodes  showed  a  noticeable  increase  in  radiation  dose  rate,  and  in  radiation  pulse  width. 
This  year  the  same  experiment  was  repeated  with  better  diagnostics,  and  with  a  heating 
power  source  that  would  work  with  one  of  the  Decade  modules,  DM2.  The  actual 
experiment  on  DM2  is  presently  scheduled  for  mid-  to  end  March. 

As  is  known,  DTRA’s  flagship  radiation  simulator  Decade  has  suffered  from  various  ills 
that  were  supposed  to  be  solved  while  the  machine  was  being  built.  However,  when  the 
pulsed  power  community  could  only  demonstrate  a  partial  solution  to  the  various 
problems,  further  construction  of  the  machine  was  suspended  until  forward  projection  of 
progress  would  give  full  confidence  that  Decade  would  achieve  its  specifications.  This 
decision  allowed  enough  time  to  take  a  new  look  at  all  the  various  components  of  a 
bremsstrahlung  machine. 

During  the  1990’s  pulsed  power  researchers  in  various  fields  had  become  convinced  that 
observed  limits  on  quantities  such  as  peak  power  and  pulse  length  in  microwave 
generators  and  the  fraction  of  lithium  atoms  in  lithium  ion  beams  were  substantially 
affected  by  surface  contamination  of  the  electrodes.  The  University  of  Michigan,  for 
example,  had  a  research  program  to  study  these  issues  in  the  context  of  microwaves.  One 
of  their  alumni,  Dr.  M.  E.  Cuneo,  managed  to  demonstrate  that  the  excess  hydrogen  in 
ion  beams  also  came  from  electrode  contamination.  This  work,  and  many  discussions 
with  Dr.  A.  Fisher,  inspired  an  SBIR  proposal  to  DTRA  by  Dr.  Pereira.  This  proposal  s 
purpose  was  to  see  whether  electrode  contamination  should  share  the  blame  for  Decade  s 
disappointing  performance  as  projected  from  experiments  on  the  two  Decade  modules 
then  existing,  DM1  and  DM2,  and  more  importantly,  to  come  up  with  a  way  to  suppress 
contamination. 

The  SBIR  Phase  I  research  successfully  demonstrated  that  heating  the  anode  to  red  and 
orange  heat,  in  the  range  of  1000  to  1200  K,  indeed  has  a  marked  effect  on  the  pulse 
length,  a  twofold  increase.  However,  the  experiments  were  done  on  a  machine,  Dr. 
Gold’s  Long  Pulse  Accelerator  at  NRL.  This  generator’s  parameters  scaled  to  DM1  and 
DM2,  but  they  were  far  from  identical.  To  make  sure  that  the  increase  in  pulse  length 
would  carry  over  to  the  Decade  program,  two  other  experiments  were  necessary.  One  was 
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a  test  of  DM2’s  plasma  opening  switch  on  DM1,  the  other  a  test  of  a  white-hot  tantalum 
anode  on  Gamble  II. 

The  main  experimenter  on  the  tests  were,  of  course,  the  scientists  who  best  know  the 
pulsed  power  machines  themselves,  Dr.  J.  Goyer  at  MPI  for  DM1,  and  Dr.  B.  Weber  at 
NRL  for  Gamble  II.  Both  tests  were  successful.  Heating  the  plasma  opening  switch 
electrodes  on  DM1  improved  the  radiation  pulse,  in  various  ways  that  are  detailed  in  a 
MPI  technical  report.  The  Gamble  II  test  developed  the  parameters  needed  for  heating  the 
tantalum  bremsstrahlung  anode  to  3000  K  without  doing  damage  to  the  remainder  of  the 
system,  and  demonstrated  also  that  the  radiation  pulse  improves.  The  results  are 
documented  in  an  NRL  technical  note,  and  will  be  written  up  in  conjunction  with  the 
experiment  that  is  now  ongoing. 

1:  The  present  effort  intends  to  test  the  white-hot  tantalum  anode  on  DM2.  The  initial 
experiment  on  Gamble  H  heated  the  tantalum  white-hot  by  a  pulse  from  a  capacitor  bank 
that  happened  to  be  available  at  NRL.  Early  on  the  best  way  to  heat  the  anode  was  not 
clear.  Continuous  heating  as  in  the  earlier  experiment  did  not  seem  particularly  practical, 
so  we  settled  on  a  heating  technique  with  an  electrical  pulse  as  first  suggested  by  NRL. 
Contrary  to  the  usual  microsecond  pulses,  the  only  pulse  that  worked  well  was  tens  of 
milliseconds  long,  at  the  limit  of  what  the  available  capacitor  bank  could  provide.  And 
the  capacitor  bank  would  be  too  small  (and  not  available)  at  DM2. 

The  proper  source  for  such  slow  pulsed  power,  as  first  suggested  by  NRL,  is  to  use  a 
series-parallel  combination  of  high  quality  lead-acid  batteries.  A  convenient  way  to 
measure  the  current  was  improvised  that  allowed  the  experiment  to  proceed  without 
hiccups.  For  the  DM2  experiment  the  current  shunt  has  been  replaced  by  a  standard  unit. 

2:  Hot  tantalum  suppresses  two  sources  of  protons.  One  is  surface  contamination,  the 
other  hydrogen  absorbed  in  tantalum’s  bulk.  Surface  contamination  can  be  handled  in 
other  ways  also,  e.  g.,  discharge  or  microwave  cleaning,  and  in  earlier  work  both  these 
techniques  were  tried.  Both  had  positive  effects.  However,  the  hydrogen  hidden  inside 
the  tantalum  stays  there.  In  fact,  tantalum  heated  to  1500  K  or  so  is  actually  a  hydrogen 
getter,  absorbing  hydrogen  from  the  vacuum  until  a  new,  lower,  pressure  equilibrium  is 
reached.  White-hot  tantalum,  with  temperature  exceeding  2300  K  or  so,  again  loses  the 
hydrogen.  The  question  is,  how  fast? 

In  the  earlier  Gamble  II  experiment  we  had  stayed  with  an  estimate  of  the  diffusion  time, 
because  the  tantalum  in  Gamble  II  is  so  thin  that  the  diffusion  time  turned  out  to  be 
almost  an  order  of  magnitude  shorter  than  the  heating  time.  However,  on  DM2  the 
tantalum  is  twice  as  thick,  0.1  mm,  and  diffusion  is  4  times  slower.  Then,  the  two 
timescales  are  more  or  less  on  the  same  order.  For  this  reason  a  more  complete 
calculation  of  hydrogen  diffusion  becomes  useful. 
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Two  independent  calculations  were  performed  giving  essentially  the  same  results.  The 
hydrogen  concentration  in  the  tantalum  drops  by  at  least  two  orders  of  magnitude  over 
the  500  ms  long  heating  pulse,  even  for  the  thicker  DM2  tantalum. 

3:  Suppression  of  the  ion  current  was  originally  measured  by  with  an  electrical  method. 
However,  the  signals  were  subject  to  interpretation.  In  an  attempt  to  diagnose  the  ions 
more  directly,  a  rudimentary  ion  mass  spectrograph  that  could  just  fit  in  the  available 
space  was  tried.  In  the  discussions  about  the  use  of  this  instrument  it  appeared  that  a 
nuclear  activation  technique  pioneered  earlier  by  NRL  would  be  easier  to  implement.  His 
nuclear  activation  measurements  demonstrate  directly  and  convincingly  that  a  white-hot 
tantalum  suppresses  the  proton  current  coming  from  the  anode,  by  more  than  two  orders 
of  magnitude  (estimated  at  400-fold). 

4;  Tantalum  stretching  on  Gamble  II  was  done  in  only  a  single  direction.  The 
experimental  procedure  is  to  clean  the  tantalum  by  a  first  heating  pulse,  and  then  shoot 
the  machine  at  the  end  of  a  second  heating  pulse.  Unidirectional  stretching  makes  the 
tantalum  wavy  in  a  single  direction,  along  the  stretch.  To  avoid  this  problem  on  DM2  a 
way  to  stretch  tantalum  in  the  radial  direction  was  designed  and  implemented.  Whether 
this  approach  is  successful  will  be  known  in  the  near  future. 
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With  sufficient  current  density,  the  electrons  in  a  bremsstrahlung  diode  can  heat  up 
the  tantalum  anode  to  such  a  high  temperature  that  the  material  evaporates  completely. 
Any  light  atoms  liberated  from  the  surface  early  in  the  pulse,  e.g.,  loosely  bound  water 
or  strongly  bound  oxides  or  nitrides,  may  then  be  swamped  by  the  hydrogen  that  was 
stored  in  the  tantalum’s  bulk.  The  surface  is  of  course  important,  and  surface  properties 
usually  dominates  in  most  common  circumstances  even  though  the  surface  contains  little 
material.  On  the  surface  there  are  light  atoms  in  perhaps  a  few  monolayers,  each  of  which 
is  maybe  0.1  nm  thick.  All  the  surface  layers  together  are  then  at  most  1  nm  or  so.  The 
tantalum  itself  is  100  000  nm  thick,  and  it  takes  only  a  very  small  fraction  of  light  ions  in 
the  bulk  tantalum  to  exceed  what  is  on  the  surface.  Therefore,  the  maximum  density  N0 
of  hydrogen  (or  other  light  atoms)  inside  the  tantalum  must  be  brought  down  to  below 
something  like  10  ppm  before  light  atoms  from  the  bulk  can  be  excluded  as  contributing 

to  a  hydrogen  (or  other  light  atom)  plasma  in  the  diode  gap. 

Earlier  experiments  have  demonstrated  that  heating  the  tantalum  close  to  its  3200  K 
melting  point  before  the  shot  has  a  strong  beneficial  effect  on  diode  behavior.  Presumably, 
the  reason  is  that  hdrogen  or  other  light  atoms  have  diffused  out  from  the  tantalum’s  bulk, 
and  also  disappeared  from  the  surface,  so  that  there  remain  only  high  atomic  number 
atoms.  These  are  heavy  and  slow,  and  it  takes  them  much  longer  than  the  pulse  time 
before  they  can  enter  the  diode  gap.  Then,  the  discharge  is  over. 

The  tantalum  in  a  continuation  experiment  is  twice  as  thick  as  the  tantalum  used 
earlier,  so  that  diffusion  is  slower  and  perhaps  too  slow.  This  note  estimates  the  diffusion 
time  in  a  one-dimensional  slab  of  tantalum  that  is  heated  by  a  slow  (30  ms)  electrica 
pulse,  in  preparation  for  the  100  000  times  faster  (100  ns)  high  power  pulse  that  follow 

after  a  few  100  ms. 

Assuming  classical  diffusion  (and  thereby  ignoring  all  kinds  of  realistic  and  very  in¬ 
teresting  and  perhaps  relevant  phenomena  such  as  the  difference  between  the  surface  and 
the  bulk  in  capturing  hydrogen,  the  effect  of  grain  boundaries,  etc),  the  one-dimensional 
linear  diffusion  equation  for  the  density  N(x ,  t)  —  Non(x,  t)  of  hydrogen  is 


^  =  D{T) *1 

dt  {  }dx2' 


(i) 


where  D(T)  is  the  linear  diffusion  coefficient  of  hydrogen  through  tantalum,  and  x  and  t 
are  the  usual  space  and  time.  The  diffusion  coefficient  D(T)  depends  on  the  temperature 
T,  but  the  temperature  is  assumed  to  be  the  same  throughout  the  tantalum  slab.  1  hen 
T  =  T(f),  and  also  the  diffusion  coefficient  is  only  a  function  of  time.  Any  skm  effect  in 
the  heating  of  tantalum  is  thereby  ignored. 
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The  normalized  hydrogen  density  n(x,  t )  vanishes  when  there  is  vacuum  on  both  tan¬ 
talum  edges,  at  d/2,  so  that  the  simplest  spatial  dependence  that  satisfies  Eq.  (1)  is 
cos(7rrr/d):  here  d  is  the  thickness  of  the  tantalum.  This  simple  cosinusoidal  profile  decays 
in  time  according  to 


dn 

dt 


D(T)tt2 


d2 


n(t), 


(2) 


where  n  =  n(t )  is  the  normalized  density.  Initially,  at  time  t  =  0,  n(t  —  0)  =  1. 

Since  the  diffusion  is  largest  at  the  highest  temperature  Tm  it  is  convenient  to  nor¬ 
malize  to  unity  according  to  D{T)  =  DmD' (T/Tm),  with  the  function  ZT(T/Tm)  =  1 
at  T  -  Tm.  However,  the  temperature  dependence2  of  the  diffusion  coeffient  follows  the 
usual  Arrhenius  factor  D(T)  =  D0  exp (~E0/kT)  up  to  the  highest  temperature  that  it  has 
been  measured  (600  K).  Extrapolating  this  relationship  to  the  much  higher  temperatures 
of  interest  here  gives 

D(T)  =  Do  exp(— Tq/T),  (3) 


where  Do/e  =  1.5  x  10-8  m2/s  is  the  diffusion  coefficient  at  the  characteristic  tempera¬ 
ture  To  =  1400  K,  so  that  Dq  =  4  x  10~8  m3/s.  Then  the  equation  for  the  normalized 
cosinusoidal  hydrogen  density  becomes 


dn 

"dt 


Dq-k2 
d 2 


exp  [— T0/T(t)]n(t). 


(3) 


Equation  (3)  is  an  ordinary  differential  equation  that  can  perhaps  be  solved  as  part  of  the 
spreadsheet  program  that  already  contains  similar  differential  equations  for  the  tempera¬ 
ture  T(f),  the  voltage  V(t),  the  current  I(t )  and  the  resistance  i2[T(t)]. 

For  estimates  to  set  the  heating  pulse  a  normalized  form  of  Eq.  (4)  is  more  appropriate, 
since  convenient  assumptions  about  the  temperature  profile  might  give  analytic  insight  in 
how  many  characteristic  diffusion  times  are  enough  to  reduce  tantalum’s  hydrogen  content 
to  an  acceptable  value.  The  normalized  form  is 

—  = -exp[-T0/T(r)]n(r),  (4) 


where  the  normalized  time  r  =  t/to  and  to  =  d2e/Don2  is  the  diffusion  time  at  the 
characteristic  temperature  To-  A  more  logical  definition  of  a  characteristic  time  would 
exclude  the  factor  e:  the  one  chosen  here  is  more  convenient,  because  the  combination 
Do/e  is  the  number  read  from  the  Arrhenius  plots. 

At  T0  =  1400  K  the  diffusion  coefficient  for  hydrogen  in  tantalum  is  Do/e  =  1.5  x  10-8 
m2/s.  For  a  0.1  mm  thick  piece  of  tantalum  heated  to  3000  K  the  diffusion  time  is  then 
t0  —  67  ms.  The  electrical  heating  pulse  in  the  previous  experiment  is  about  30  ms,  but 
for  that  experiment  the  tantalum  was  half  as  thick  and  the  diffusion  time  was  17  ms,  four 
times  shorter  than  now.  Certainly  then,  but  also  now,  the  diffusion  time  is  short  enough 
for  the  density  to  decay  appreciably,  especially  since  the  temperature  decays  slowly  after 
heating.  The  total  time  that  the  tantalum  is  hotter  than  T0  =  1400  K  may  be  a  few  100 
ms,  so  that  it  seems  safe  to  ignore  the  details  as  expressed  by  Eq.  (4). 
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